^ 5?V, 3-2-5 

(12) INTERNATIONAL yi^^ICATION PUBLISHED UNDER THE PATENT (^S'ERATION TREATY (PCT) 

IlllliiOlllllllllllllllililllllllllllllllllili 

(10) Internatioiial Publication Number 

wo 2004/015631 A2 



(19) World Intellectual Property 
Organization 
International Bureau 

(43) International Publication Date 
19 February 2004 (19.02.2004) 




PCT 



(51) International Patent Classification'': 



G06T 11/00 



(21) International Application Number: 

PCT/CA2003/001214 

(22) International Filing Date: 8 August 2003 (08.08.2003) 

(25) Filing Language: English 

(26) Publication Language: English 



(30) Priority Data: 
2.397,389 



9 August 2002 (09-08.2002) CA 



(71) Applicants (for all designated States except US): UNI- 
VERSITE DE SHERBROOKE [CA/CA]; 2500. boule- 
vard de rUniversite. Sheibrooke. Quebec JIK 2R1 (CA). 
BISHOP'S UNIVERSITY [CA/CA]; Lennoxville, Que- 
bec JIM 1Z7 (CA). 

(72) Inventors; and 

(7^ Inventors/Applicants (for US only): ZIOU, Djemel 



[CA/CA]; 2665 Maricourt, Sherbiooke, Quebec JIK 1R6 
(CA). AUCLAIR-FORTIER, Marie-Flavle [CA/CA]; 
264 Lockwell, Quebec, Quebec GIR 1 V9 (CA). ALLILI, 
Madjid [CA/CA]; 585 London, apt. 1. Sherbrooke, 

Quebec J IH 3N2 (CA). 

(74) Agents: BROUILLETTE, Robert et al.; Brouillette 
Kosie Prince, 1100 Ren6-L^vesque Blvd. West. 25th 
Floor, Montreal, Quebec, H3B 5C9 (CA). 

(81) Designated States (national): AE, AG, AL, AM, AT, AU, 
AZ, BA, BB, BG, BR, BY. BZ, CA, CH, CN, CO, CR, CU, 
CZ. DE. DK, DM, DZ. EC, EE. ES. FI, GB, GD, GE. GH, 
GM, HR, HU. ID. IL. IN. IS, JP. KE, KG, KP, KR, KZ, LC. 
LK, LR, LS. LT. LU, LV, MA. MD, MG, MK, MN, MW, 
MX, MZ, NI, NO, NZ, OM, PG, PH, PL, PT, RO, RU, SC, 
SD, SE, SG, SK, SL, SY, TJ, TM, TN, TR. TT, TZ. UA, 
UG, US, UZ, VC, VN, YU, ZA, ZM, ZW. 

(84) Designated States (regional): ARIPO patent (GH, GM, 
KE, LS, MW, MZ, SD, SL. SZ, TZ, UG, ZM, ZW), 

[Continued on next page] 



(54) Title: IMAGE MODEL BASED ON N-PIXELS AND DEFINED IN ALGEBRAIC TOPOLOGY, AND APPLICATIONS 
THEREOF 



DEFINING A PIXEL DIMENSION n TO PRODUCE 
A n-PIXEL IMAGE SUPPORT 




r 


EXPRESSING THE rvPlXELS OF THE IMAGE 
SUPPORT ALGEBRAICALLY IN RELATION TO 
THE ivPIXEL DIMENSIONS 


1 


r 



PRODUCING A GEOMETRICAL COMPLEX 
INCORPORATING AT LEAST ONE ivPIXEL 
IMAGE SUPPORT 



6601 



6602 



6603 



(57) Abstract: A computational image model comprises an image 
support including a structure of n-pixels comprising pixel faces, 
quantities related to image features, and an algebraic structure 
relating the quantities to the n-pixels and/or pixel faces, the 
algebraic structure comprising algebraic operations defining a 
relation between the quantities. A method of computationally 
modelling an image comprises producing an image support 
including a structure of n-pixels comprising pixel faces, defining 
quantities related to image features, and relating the quantities to 
the n-pixels and/or pixel faces through an algebraic structure, and 
relating the quantities to each other through algebraic operations. 



o 



EXPRESSING THE GEOMETRICAL COMPLEX 
ALGEBRAICALLY BY MEANS OF q<^HAINS 



EXPRESSING THE SCALAR VECTORS AND 
MATRICES BY MEANS OF THE COEFFICIENTS 
OF THE q^HAlN 



6604 



6605 



EXPRESSING GLOBAL QUANTITIES WITH ALL 
q-PIXELS THROUGH q-COCHAINS Fq 
COMPRISING ASSOCIATING GLOBAL 
QUANTITIES TO THE Q-PIXELS AND ITS FACES 



6606 




104/015631 A2 



Eurasian patent (AM, AZ, BY, KG, KZ, MD, RU, TJ, TM), 
European patent (AT, BE, BG, CH, CY, CZ, DE, DK, EE, 
ES, H, FR, GB, GR, HU, IE, IT, LU, MC, NL, PT, RO, 
SE, SI, SK, TR), OAPI patent (BF, BJ, CF, CG, CI, CM, 
GA, GN, GQ, GW, ML, MR, NE, SN, TD, TG). 

PubUshed: 

— without international search report and to be republished 
upon receipt of that report 




For two-letter codes and other abbreviations, refer to the "Guid- 
ance Notes on Codes and Abbreviations" appearing at the begin- 
ning of each regular issue of the PCT Gazette. 



^|P^CT/CA2003/001214 

, :^(^ -;FTO ^FEB 200S 



IMAGE MODEL BASED ON N-PIXELS AND DEFINED IN 
ALGEBRAIC TOPOLOGY, AND APPLICATIONS THEREOF 



FIELD OF THE INVENTION 

The present Invention relates to an image model based on n-pixels. 
More specifically, the present Invention is concerned with an image model 
10 based on n-pixels, defined in algebraic form and applicable, in particular but 
not exclusively, for the resolution of diffusion and optical flow, and for the 
deformation of curves. 



15 BACKGROUND OF THE INVENTION 

People have a notion of what an image Is. For instance, for 
psychologists the Image is linked to the shape of objects, their depth, the 
relationship of these shapes and their piBrceptual organization. 

20 

Artists are focused on how features such as shape, color, and 
perspective are organized to represent a scene that may originate in their 
imagination. 

25 Physicists are concerned with the physical phenomena produced by a 

given scene and how they are represented in the image. 



30 



Neurophysiologists regard images through visual phenomena in 
humans and animals, such as contrast sensitivity, Mach bands, contrast, 
constancy, and depth perception. 
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While a precise definition of the image is elusive, it is clear that for 
certain people an image is the visualization of physiological, perceptual, or 
physical phenomena and for others it is related to semantic content. Whatever 
the term image means, the intention Is to establish a foundation upon which 
5 images of all forms and contents can be discussed with minimal confusion. 

In image processing and computer vision, the foundation is related to 
the understanding of the image formation process. Light generated by a 
source is modified by the objects of the scene. The modified light is captured 

10 by a system acquisition device, transformed into an appropriate form and 
displayed on a physical support. The content of the resulting image and, 
consequently, its further use, both depend on the properties of the light 
(structured or not, spectral properties such as the range of wavelengths, the 
number of ranges and the intensity), the characteristics of the acquisition 

15 device, the transformation (analog-to-digital, pre-processing), the display 
format (temporal or spatial organization of image elements, film, vector, 
raster). From the automatic processing point of view, the image is enhanced to 
improve its perceptual quality or to make some of its features explicit. Usually, 
its content is analyzed via a successive reduction process to construct a more 

20 descriptive representation in terms of relevant features, which can be used 
more effectively by a decision system (car, robot, etc.) or to help human 
beings in their daily tasks. 

One of the concerns here is to focus on the data structure for images 
25 and its consequences on the processing scheme. 

Algebraic topology concepts are a key to representing images. 
Algebraic topology is well-known domain in mathematics [10, 6, 9]. The 
literature of algebraic topology offers wide knowledge that can be used for 
30 Images. The specialists for algebraic topology, however, have made no effort 
to implement their knowledge into computer vision and image processing. 
Instead of using algebraic topology, the specialists of image processing and 
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computer vision have limited themselves to develop solutions based on the 
topology and on discrete geometry [7, 8], 

Algebraic topology was first introduced Into image processing and 
5 computer vision by Allill and Ziou for topological feature extraction and shape 
representation In binary Images [1, 2, 15]. This technology is used by Allill, 
MIschaikow and Tannenbaum [3] in medical binary images. Audafr et al. [4] 
also used algebraic topology for linear and non-linear isotropic diffusion as well 
as for optical flow in gray level images. P. Poulln et al. [12] used algebraic 
10 topology for snakes and elastic matching In gray level images. 

In image processing and computer vision, several Image models have 
been accepted and are in recurrent use since several decades. These image 
models integrate both the image support and the local quantities associated 

15 with this support. The Image support Is fonned by pixels. With each pixel. Is 
associated a scalar quantity called a gray level, or a vector quantity called 
either color at the perceptual level or multispectral at the signal level. Existing 
models differ primarily in terms of how the image support (definition and the 
organization of pixels) and how values are fomnulated. The well-known image 

20 models are the function, the random process, and the ordered set. The image 
is a function x -» G"' , where L, = {l,...,N^} and = {!,..., N^} , N.xN, is 
the resolution of the image, G = {1, ... ,n}, where n Is the maximal quantity and 
m the number of image bands. In the case of a binary Image n = 2, image 
processing has roots In graph theory, language theory, logic, discrete 

25 geometry, and so on. If n > 1, usually the image is modelled as a real function 
(analogue image where G,L^,Ly <z 91 ). In this case, function theory, functional 
analysis, catastrophe theory, differential equations, and differential geometry 
constitute the foundation. An image can also be modelled as a collection of 
random variables -?"(/, y)|(z,y)eX,xX^. In this case, the probability density 

30 function, moments, sufficient statistics, time series, and the Markov processes 
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are the roots. When the image Is modelled as an ordered set, discrete 
mathematics and mathematical morphology are the foundation , 

Since the introduction of mathematical morphology, efforts of 
5 researchers in these fields have been focused on the use of more and more 
complex mathematical, physical or computer concepts as fonmalism of specific 
problems (edge detection, Image segmentation, optical flow, and defomiation) 
without questioning the image model. The definition of a new image model can 
lead algorithms that are designed, on new basis. An image is a physical or 
10 mathematical quantity where variables (image support) represent geometrical 
or temporal elements such as points, lines, surfaces, and times. For example, 
the image, as a function xZ^ -> G"* , can be defined by both the geometrical 
and topological properties of the domains L^xLy and the topological, 

geometrical and analytical properties of C?". Although existing image models 
15 have deep roots in mathematics or In physics, the variables, the quantity and 
the association between them are not well-defined. For a given computer 
vision or Image processing task, no formal mechanism is given for the 
integration of physical, topological, geometrical properties of objects and their 
behaviours as a part of the image model. Consequently, the resulting 
20 computational schemes are non-modular and sometimes not easy to 
reproduce. 

SUMIVIARY OF THE INVENTION 

25 

According to the present invention, there is provided a computational 
image model, comprising an image support including a structure of n-pixeis 
comprising pixel faces, quantities related to image features, and an algebraic 
structure relating the quantities to the n-plxels and/or pixel faces, the algebraic 
30 structure comprising algebraic operations defining a relation between the 
quantities. 



wo 2004/015631 




•CT/CA2003/001214 



The present invention also relates to a method of computationally 
modelling an image, comprising producing an image support including a 
structure of n-pixeis comprising pixel faces, defining quantities related to 
5 * innage features, and relating the quantities to the n-pixels and/or pixel faces 
through an algebraic structure, and relating the quantities to each other 
through algebraic operations. 

Still in accordance with the present invention, there is provided a 
10 computational framework for solving a heat transfer problem, comprising: . 

producing an image support including a structure of n*pixels, the image 
support comprising: 

q-pixels respectively translating the n-pixel algebraically, wherein q 
e {1, 2,..., n}, and wherein each q-plxel includes (q-l)-faces, (q-2)- 
15 faces, (q-q)-faces; 

geometrical complexes each being a collection of q-pixels; 
q-chains respectively expressing the geometrical complexes in 
algebraic form, each q-chain being a linear combination of all the q- 
pixels of the geometrical complex; 
20 - in the geometrical complexes, q-cochains which are relations 

associating quantities related to image features to the q-pixels 
and/or faces of said q-pixels; and 
a coboundary defining a relation between q-cochains; 
computing a q-cochain T of a first geometrical complex as the location 
25 of unknown temperatures; 

corinputing a q-cochain H of the first geometrical complex as a global 
temperature variation; 

finding a q-cochain 8 of a second geometrical complex as a global 
energy variation, as a function of the q-cochain H through a linear 
30 transformation; 

finding the q-cochain e as a function of the q-cochain T; 
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defining a q-cochain G of the first geonnetrical complex from the q- 
cochain T through a first coboundary operation, transforming the q-cochain G 
into a q-cochain Q of the second geometrical complex, and defining, from the 
q-cochain Q and through a second coboundary operation, a q-cochain D of the 
5 second geometrical complex as a global diffusion; 

defining a q-cochaIn S of the second geometrical complex as a global 
source; 

establishing a relation between the q-cochains s, D and S. 

10 The present invention further relates to a computational framework for 

two-dimensional active contour model, comprising: 

producing an image support Including a structure of n-pixels, the image 
support comprising: 

- q-pixels respectively translating the n-pixel algebraically, wherein q 
15 € {1, 2,..., n}, and wherein each q-pixel includes (q-l)-faces, (q-2)- 

faces, (q-q)-faces; 

- geometrical complexes each being a collection of q-plxels; 

- q-chains respectively expressing the geometrical complexes in 
algebraic form, each q-chain being a linear combination of all the q- 

20 pixels of the geometrical complex; 

in the geometrical complexes, q-cochains which are relations 
associating quantities related to image features to the q-pixels 
and/or faces of said q-pixels; and 

- a coboundary defining a relation between q-cochains; 

25 computing a displacement q-cochain D of a first geometrical complex; 

computing a strain q-cochain S of a second geometrical complex, 
comprising: 

- defining an approximate strain function e(x) as a function of 
the q-cochaIn D; 

30 - expressing the q-cochain S as a function of the approximate 

strain function and relative positions of the first and second 
geometrical complexes; and 
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computing a force q-cochain F of the second geometrical complex as a 
coboundary of the strain q-cochain S. 

The foregoing and other objects, advantages and features of the 
5 present invention will become more apparent upon reading of the following 
non-restrictive description of illustrative embodiments thereof, given by way of 
example only with reference to the accompanying drawings. 

The present specification refers to various references. These 
10 references are herein incorporated by reference in their entirety. 

BRIEF DESCRIPTION OF THE DRAWINGS 

15 In the appended drawings: 

Figure 1 Is an example of a 2-cube; the orientation is given by definition 
[0.1] X [0,1]: 

20 Figure 2a is an example of subdivision that is a cubical complex, and 

Figure 2b is another example of subdivision that is not a cubical complex; 

Figure 3 illustrates, in solid lines, an example of a primary cubical 
complex and, in dashed lines, an example of a secondary cubical complex; 

25 

Figure 4a is an example of original image and Figure 4b is an example 
of a resulting smoothed image; 

Figures 5a, 5b and 5c illustrate a body in 3D space at time t. In Figure 
30 5a, xi is the location of the particle Xi, n is a vector normal to the surface, dS 
is an infinitesimal surface patch and dV is an infinitesimal amount of volume. In 
. Figure 5b, fe is an external force applied on dS, p is the mass and b is a body 
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force applied on dV. In Figure 5c, q is the heat flow passing through dS and r 
is the heat produced by dV; 

Figure 6a, 6b and 6c are three examples of q-pixels in (h x I2). 
5 Figure 6a illustrates an example of 0-pixel where I1 = {2} and I2 = {1}, Figure 6b 
an example of 1 -pixel where I1 = [2,3] and I2 = {1}, and Figure 6c an example 
of 2-pixel where I1 = [2,3] and I2 = [1,2]; 

Figure 7 is an example of coboundary operation; 

10 

Figure 8a is an example of projection orito the tangential part of the 
domain, and Figure 8b is an example of projection onto the normal part of the 
dorhain; 
• 

15 Figures 9a and 9b are examples of cochains for one 3-pixel of 

Figure 9a illustrates an example of cochain T while Figure 9b illustrates an 
example of cochain H; 

Figure 10a is an example of cochain Q for one 3-pixel of K®, and Figure 
20 1 0b is an example of cochain D for one 3-pixel of K^• 

Figure 11a is an example of cochain T for one 3-pixei of K*', and Figure 
1 1 b is an example of cochain G for one 3-pixel of K^; 

25 - Figure 12 is an example of three different paths between two points; 

Figure 13 is an example of computational scheme for an unsteady 
problem with no source; 

SO Figure 14 is an example of two cubical complexes for a 5 X 5 Image; 

Figure 15 is an example of yp for one 2-pixel of K^; 



wo 2004/015631 




►PCT/CA2003/001214 



9 



Figure 16a Is an example of ye (in dashed lines) for one 2-plxel of 
Intersecting four pixels of K^, Figure 16b is an example of cochains T and G, 
and Figure 16c is an example of cochain Q; 

5 

Figure 17 is an example of one 3-pixel of IC' surrounding one 1-pixel of 
Figure 18 is an example of Xs for one 2-pixel of K^; 

10 

Figure 19a, 19b and 19g are examples of one 3-pixel of K® intersecting 
four 3-pixels of K^, for cochain T (Figure 19a), cochain G (Figure 19b), and 
cochain Q (Figure 19c); 

1 5 Figure 20 is an example of two 2-pixels of with Ts; 

Figure 21a Is an example of physics-based isotropic diffusion o = {2.0, 
4.0, 5.0}, and Figure 21b is the example of physics-based isotropic diffusion of 
Figure 21a convolved, for same scales; 

20 

Figures 22a, 22b and 22c are examples of first images of three 
sequences, a rotating sphere sequence (Figure 22a) a Hamburg taxi 
sequence (Figure 22b) and a tree sequence (Figure 22c); 

25 Figure 23a is an example of flow pattern computed for the rotating 

sphere sequence of Figure 22a using a physics-based nnethod, and Figure 
23b is an example of flow pattern computed for the rotating sphere sequence 
of Figure 22a using an iterative finite difference method; 



30 Figure 24a is an example of flow pattern computed for the Hamburg 

taxi sequence of Figure 22b using the physics-based method, and Figure 24b 
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is an example of flow pattern computed for the Hannburg taxi sequence of 
Figure 22b using the iterative finite difference method; 

Figure 25a is an example of flow pattern computed for the tree 
5 sequence of Figure 22c using the physics-based method, and Figure 25b Is an 
example of flow pattern computed for the tree sequence of Figure 22c u^ing 
the iterative finite difference method; 

Figure 26a is an example of a first image of the tree sequence of Figure 
1 0 22c with white noise added (standard deviation of 1 0); 

Figure 27a Is an example of flow pattem computed for the tree 
sequence of Figure 22c with white noise added using the physics-based 
method, and Figure 27 b is an example of flow pattem computed for the tree 
5 sequence of Figure 22c with white noise added using the iterative finite 
difference method; 

Figure 28a Is a secUoh of the peppers Image (o = 5) of Figure 21a 
corresponding to the original Image with noise added, Figure 28b is a section 
20 of the peppers image (a = 5) of Figure 21a obtained with the PB method, and 
Figure 28c is a section of tlie peppers image (o = 5) of Figure 21a obtained 
with the FD method; 

•Figure 29a, 29b, 29c and 29d are examples of a section of the peppers 
25 image (o = 5) of Figure 21a corresponding to the original image (Figure 29a), 
corresponding to the original with white noise added (Figure 29b, obtained 
using the PB method (Figure 29c), and obtained with the FD metliod (Figure 
29d); 

30 Figure 30a is an example of PB method for o = {1 .0, 3.0, 5.0}, and 

Figure 30b is an example of FD method for the same scales; 
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Figure 31a is an example of original image, and Figure 31b . is an 
exanriple of image witli noise added (standard deviation =10); 

Figure 32a is an example of PB method for a = {4.0, 8.0}, and Figure 
5 32b is an example of FD method for the same scales; 

Figure 33a is an example of PB method, and Figure 33b is an example 
of FD method with o = 8.0; 

10 Figure 34 is an example of a body of arbitrary size, shape and nnaterial, 

where AS is a surface patch, f is a vector of surface forces applied on AS, AV 
is an amount of volume with a mass p, and b is a vector of body forces applied 
on AV; 

15 Figure 35 is an example of cutting plane passing through a point and 

partitioning the body into two sections; 

Figure 36 is an example of a force Af acting on a cutting plane with 
normal vector n; 

Figures 37a is an example of stresses on the original body, Figure 37b 
is an example of normal stress after an extension of the body, Figure 3Tc is an 
example of normal stress after a compression of the body, and Figure 37d is 
an example of shear stress after a distortion of the body; 

Figure 38 is an example of the component of the stress in the direction 

of xi; 

Figures 39a and 39b are examples of the deformation (Figure 39a) and 
30 distortion (Figure 39b) of a body subjected to stresses; in both Figures 39a 
and 39b, the rectangle ABCD has been deformed or sheared into A'B'CD; 



2b 



25 
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Figure 40 Is an example of normal strain of a body; 
Figure 41 is an example of shear strain in a body; 

Figure 42 is an example of a body B in motion and subjected to forces, 
5 wherein ti" and pbj are respectively the traction and body forces in the direction 

OfXii 

Figure 43a is an example of kinematic equation. Figure 43b is an 
example of constitutive equation, and Figure 43c is an example of 
10 consen^ation equation; 

Figure 44 is an example of decomposition of the linear elasticity 
problem into basic laws; 

15 Figures 45a, 45b and 45c are three examples of q-pixels in (li x b). 

More specifically, Figure 45a is an example of 0-pixeI for h = {2} and I2 = {1}, 
Figure 45b is an example of 1 -pixel for U = [2,3] and I2 = {1}, and Figure 45c is 
an example of 2-pixel for h = [2,3] and I2 = [1 .2]; 

e 

20 Figure 46 is an example of the coboundary operation; 

Figure 47a is an example of cochain u for a 3-plxel of K^. and Figure 
47b is an example of cochain D for a 3-pixel of K^; 

25 Figures 48a and 48b are example of cochains S and F, respectively, for 

a 3-pixel of K^; 

Figure 49 is an example of two cubical complexes for a 5 x 5 image; 

30 Figure 50 is an example of a 2-pixel of and the topological quantities 

associated with it; 
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Figure 51a is an example of yf in dashed lines. Figure 51b is an 
example of 2-cochains D and u, and Figure 51c is an example of 2-cochain S; 

5 Figure 52 is an example of five adjacent vertices in a curve and its 

deformation; 

Figure 53 is a table that shows typical values 6f the Poissoh's ratio and 
Young's modulus of elasticity for some materials; 

10 

Figure 54a is an example of initial curve, and Figure 54b is a bright line 
plausibility image; 

Figure 55a is an example of results obtained for an aerial image using 
15 the PB method, and Figure 55b is an example of results obtained for an aerial 
image using the FEM method; 

Figure 56 is an example of initial curve for a SAR image; 

20 Figure 57 is an example of line plausibility for a SAR image; 

Figure 58 is an example of road correction for a SAR image with the PB 
method; 

25 Figure 59 is an example of road correction, for a SAR image with the 

FEM method; 

Figure 60 is an example of initial curve; 



30 



Figure 61 is an example of bright line plausibility image; 
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Figure 62 is an example of corrections for a multiband Image (PB 
metliod); 

Figure 63 is an example of corrections for a Landsat 7 image (FEIVI); 

5 

Figure 64a is an exampie of initial curves for a synthetic image, and 
Figure 64b is an example of corrected curves for a synthetic image; 

Figures 65a, 65b and 65c show an example of shape recovery of a 
10 curve when the external forces are removed, and Figure 65d is an example of 
both Initial curve (in white) and final curve (in black); 

Figure 66 is a schematic flow chart showing how to build an illustrative 
embodiment of the computational image model according to the invention; and 

Figure 67 is a schematic flow chart showing how to biiild a 
computational framework for solving a problem, in accordance with an 
illustrative embodiment of the present invention and using the computational 
image model of Figure 66. 

DETAILED DESCRIPTION OF THE ILLUSTI^TIVE EMBODIMENTS 

The illustrative embodiments of the present invention are generally 
25 based on a data structure for images based on n-pixels, in which the image 
support, the image quantities and the allowable operations are specified 
separately. In this data structure, mathematics and physics are unified; that is, 
the data structure allows taking into account constraints originating in physics, 
mathematics, and the future use of the Image. The image dimension is explicit, 
30 which allows us to design algorithms that operate in any dimension. The 
foregoing and other advantages and new open problems of this data structure 
will be discussed herein below. 



15 



20 
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One of the goals of the present invention is to provide a computational 
image model or a data structure that is capable of capturing all object 
properties that are. needed for a given task. Figure 66 is a schematic flow chart 
5 summarizing the procedure (successive steps 6601-6606) for building such a 
computational image model according to an illustrative embodiment of the 
invention. A data structure of the image is the fomial specification of the image 
variables, image quantities, the association between image quantities and 
variables that enables capture of the geometrical and topological properties of 

10 objects as well as their physical and mathematical behavior. The abstract view 
of an image as a data structure as is used in computer programs is defined by 
the attributes and a collection of meaningful operations. The attributes are the 
image, support and quantities that are assigned to the image support such as 
the image radiometry (e.g., color and grey level) or any feature that can be 

15 deduced from the radiometry (e.g. texture). These quantities are scalar, vector 
or tensor. The allowable operations are of two kinds: the operations that are 
problem independent such as read and write and those that are problem 
dependant such as object deformation, diffusion and optical flow. To 
summarize, the image is defined by the support, quantities and allowable 

20 operations. The iatter three iterhs will be defined in detail in the following 
description. 

Image support 

25 Often, discretization of an image is achieved via a cubic tessellation. 

For example, a 2D (two-dimensional) image support is formed by unit squares 
commonly called pixels. Similarly, a 3D (three-dimensional) image is a 
tessellation of unit cubes commonly called voxels. More generally, the 
illustrative embodiments of the present invention will consider the image in n 

30 dimensions as a set of unit n-cubes, which are commonly called n-pixels. 
When n = 0, the image is a set of points; when n = 1 , a set of edges; when n = 
2, a set of squares; when n = 3 a set of cubes, and so on. Any two n-pixels are 
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either disjoint or intersect througli a comnnon i-pixel, where i < n. This 
subdivision of the image support is not unique. Several other geometrical 
forms such as, for example, triangles or hexagons can be used. It has been 
proven that the topological features of the image support do not depend on the 
5 subdivision used [13]. The cubical subdivision is commonly used in image 
processing and computer vision and will therefore be used in the non- 
restrictive illustrative embodiments of the present invention. One does hot 
need to explicitly define an orientation of pixels. Indeed, the definition of a pixel 
as a product of intervals in a certain order coupled with the natural order of 
10 real numbers imposes an orientation on each coordinate axis and also on the 
canonical basis of SR" (see Figure 1 ). 

Formally, a pixel ctg 91" as a geometric entity is translated into an 
algebraic structure as follows: 

15 

. " (1) 

where x is the Cartesian product and li is either a singleton or an interval of 
unit length with integer endpoints; i.e., 1] is either the singleton {k} , in which 
20 case it is said to be a degenerate interval, or the closed interval [A:,Ar + l] for 
some integer k. The number ^ e {O,l,...,n}of non-degenerate intervals In 
Equation (1) is, by definition, the dimension of cr , which will be referred to as a 
q-pixel. For q>l, let 7 = {ko,ki,...,kq.i}be the ordered subset {l,2,...,n}of 

indices for which 4^ = l^jybj] is a non-degenerate interval. Define 



25 



Ak^a = /i X ... X Ikj^i X {aj} x 4^.4.1 x . . . x 4 



and 



wo 2004/015631 



17 



IPCT/CA2003/001214 



Bk^a = /i X • . . X 4^„i X {6j} X /fc^.+i X . _ X 4 

where ^^^cr and J?jt;Cr are, respectively, the front (q-l)-face and the back (q-1)- 

face of o-. Each of these (q-l)-faces is a (q-l)-pixels. These faces are then 
5 called (q-l)-faces of cr . In the sanrie manner, one can define the (q-2)-faces, 
, down to the 0-faces of a . Figure 1 shows a 2-pixei A, with its 1 -faces a, 
b, c, d. The 0-faces of A are the vertices that are not represented here for the 
sake of clarity of the picture. The boundary of the q-pixels a enables the 
writing of the relationship between a q-plxels and Its (q-l)-faces in algebraic 
10 fonri. By definition this is the alternating sunn of its (q-1 >faces, i.e. 

fco (2) 
For example, the boundary of A in Figure 1 1s given by the following 

15 relation: 

d^a = (-0 X [0,1] + 1 X [0,1]) - ([0,1] X 0 - [0,1] x 1) = (c - a) + (c/ ^ 6) . 

So far, the pixel, its faces, and the association between them have 
20 been defined geometrically and algebraically. Now, tine image support will be 
defined as a geometrical entity, called a cubical complex, and then its 
algebraic structure will be written. A cubical complex K in 51" is a collection of 
q-pixels where 0^q<n such that : 

25 - Every face of a pixel in K is also located in K; and 

- The intersection of any pair of two pixels of K is either empty or formed 
by a common face of both pbcels of the pair. 
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The first condition implies that all faces of a pixel belong to the cubical 
complex. 

The second condition concems the organization of the cubical complex. 

5 If the intersection of two pixels is empty, then the image support is formed by 
one or more connected components. This condition provides an image support 
that is more general than existing definitions since it allows a formal 
specification of a cubical complex that is formed either by one or more image 
supports (e.g., an image sequence) or by several distinct binary objects. When 

10 the intersection is a face, certain geometrical configurations of the complex are 
ruled out. For example, Figure 2a illustrate a subdivision that is a cubical 
complex and Figure 2b illustrates a subdivision that is not a cubical complex. 

The dimension of the cubical complex K is, by definition, the largest 
15 number q for which K contains a q-pixej. 

As in the case of the q-pixel, the cubical complex can be written In 
algebraic form. Given a topological space X <^^R" in terms of a cubical 
complex, the set of all q-pixels of X is denoted ={o-q,.,.,crJ^^}. A.q-chain In 
20 X is a formal sum of integer multiples of elements of . More precisely, It is a 
linear combination 

25 where Ai,...,A^^are integers. For example, in Figure 1, a-c + d-b is a 1- 
chain. Two q-chains are added by adding corresponding coefficients. The set 
of q-chains can be given the structure of a free Abelian group with basis , 

usually denoted by C^iX), Since only finite complexes will be considered, the 



wo 2004/015631 ^j^j^, ^^CT/CA2003/001214 

19 



groups C^(X) are finitely generated and C^(Ar) = 0 if q is greater tlian the 
dimension of the complex; naturaily, C^iX) = 0 If q < 0. 



To define the chains that are associated with tlie faces of pixels of a 
5 cubical complex. Equation (2) is extended by linearity to all q-chains to obtain 
the boundary map d^ :C^(X)^C^_^(X) . Note that do = 0 since C^i(A^) = 0- 
The boundary map satisfies the following property [9]: 

dq o dq+i = 0 

10 

To summarize, the discrete image support of any dimension n is 
formed by n-pixels. Unlike conventional image models, the n-pixel is . a 
dimensional geometric entity formed by other geometric entities called faces, 
The geometrical n-pixel is translated into a recurrent algebraic structure, more 

15 reliable for mathematical handling. All n-pixels of the image support form a 
cubical complex, a geometrical entity that is translated into an algebraic 
structure called the chain. The association between the n-pixel and its faces, 
and hence between chains of successive dimensions is given by a boundary 
operator. It should be noted that the use of any other geometrical primitive 

20 such as a triangle or hexagon for subdividing the image support affects the 
computational complexity of the derived algorithms, but it affects neither the 
topological features of the image support nor the .computation rules for these 
features. 

25 Image quantities 

In the foregoing description, the concept of the finite cubical complex 
was introduced to give an algebraic description of the discrete image support. 
A similar formulation is needed to describe the image field (scalar, vector, 
30 matrix) over the discrete image support. For this purpose, let's return to the 
above described chain concept to give a more general definition. Considering 
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a cubical complex K of dimension n, each q-plxel (9 ^ 71 ) of K is associated to 
a coefficient in the ring (A,+ *), where the elements may be scalars (gray 
level), vectors (color, gradient), matrices (Hessian), etc. The chain is the 
formal sum ^^,4 ' where ;i. is a coefficient in (A,+,*), and the generators c\ , 
5 Vz form a basis of an Abelian group. The chain can be seen as a vector 
where iST^is the number of generators used. Consequently, two 

chains can be summed by adding their coefficients and multiplied using the 
scalar product. The addition and multiplication are taken according to the 
definition of a ring. Moreover, one can define a null chain (respectively unit 
1 0 chain) whose coefficients are all equal to the null (respectively unit) element of 
the + (respectively *) operation of (A,+,*). Consequently, q-chains define an 
image model that has attractive computational properties since they form a 
rich algebraic structure, the module (i.e, a vector space defined on a ring). 



To briefly show how to use the chain nnodel in image processing, a 
simple illustrative example concerning any global transform of the image such 
as histogram equalization or thresholding is presented. The ghain coefficients 
are the gray levels. The formal expression' of the global transform implies two 

applications: H : (^^,+,*) and H : (^,+,*) (^,-1-,*) , where (^^,+,*) 

is a module. They are defined by Hi^XiC^ = X/^^'^/)^*' • 

The drawback of chain-based image models is that the physical or 
mathematical quantities and the image support are described together in a 
formal sum. Consequently, the chain coefficients combine mathematical or 
25 physical quantities, pixel orientation, and possibly other quantities such as 
weights associated with pixels. For example, there is ambiguity in the 
interpretation of the sign of it may be due to the orientation, the physical 
quantities, the weights or their multiplication. This image model can be 
considered adequate, especially in physics [14], engineering and computer 
30 graphics [1 1 , 5], where chains have been used to model quantities. However, 
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this illustrative embodiment of the present invention proposes to refine this 
quantity model to overcome the confusion produced by the chain coefficients. 

To reflect only the geometry of the image support (e.g., orientation and 
5 multiplicity), what follows will consider q-chains with integer coefficients. We 
are looking for an application :C^(^-->(jB,+,*), which associates a global 

quantity (energy, gray level, color, flux, tensor, etc.) with all q-p!xels, where 
q<.n and (5,+,*) is a ring. As in the case of the chain-based image model, for 

two adjacent q-pixels c\ and c^, must satisfy 

10 +>^7C^) = >^^,(«?J) + >l2^^(c^), which means that the sum of the 

quantities generated within each q-pixel is equal to the quantities generated 
within the two q-pixels. F^ can be extended by linearity to any q-chain ^^fi^ , 

where each is an integer, as follows: 

F^\s called a q-cochain. To illustrate the cochain concept, let us 
consider a 2-pixel oz and a vector field V. A 0-cochain is defined by the value 
of V at 0-pixels. A 1 -cochain is i the line integral over the faces of ca- A 

20 2-cochain is ^ dS , the surface integral over the 2-pixel C2, and the the dot 
product. 

To summarize, a cochain allows us to associate quantities with the q- 
pixels and with the faces thereof. Unlike existing image models, the model 
25 according to the illustrative embodiments of the present invention provides a 
rich structure that allows the definition of both local and global quantities. 
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Operations 

A generic operation that can be instantiated depending on the problem that 
' we are dealing with will now be . defined. Having in mind that a q-pixel has 3^ 
5 faces, the generic operation should specify the algebraic relationship between the 
quantities (i.e., cochains) associated with these faces. Based on the linearity 
principle, the quantity of a given q-pixel is transferred to its cofaces with the same 
or opposite sign, according to the agreement of its orientation and the orientation 
of its cofaces. The quantities that are transferred to the (q+1)-pixels by its faces 
10 are summed. More formally, it should be reminded that the relationship between 
the (q-l)-chain and the q-chain is given by the boundary operator. Similarly, the 
relationship between the q-cochain and the (q+1)-cochaln is given by the 

coboundary operator S^:C^ C^*^ , where C is the Abelian group of q- 

cochains. Given a (q+1)-chain c, this coboundary operator Is defined by: 



The capacity of the cochain and the coboundary to model a given problem 
will now be discussed. The cochain is a linear application and should fulfill the 

20 coboundary requirement of Equation (5). Thus a question concerns the limits 
in modeling a given quantity. It is difficult to answer this question for the 
general case. Much investigation is needed first. In the present illustrative 
embodiment, this question is only answered by identifying several problems 
that can be modeled by the cochain and coboundary. Certain problems such 

25 as convolution and its applications (smoothing, numerical differentiation, high- 
pass filtering, noise estimation) and the Fourier transform can be modeled by 
the cochain without approximation since they only require setting the 
coefficients /L^ of the cochain to specific values (see [16] for the case of 
numerical differentiation). Thresholding can be represented by the cochain 

30 without approximation since ^(^/l-ifi) = » 

whereiir:(5,+,*) ->(5,+,*) . Other problems can be broken down to basic 



15 




(5) 
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laws, each of which can be. described by the topological Equation (5). For 
example, it has been pointed out [14] that many physics problems can be 
broken down into basic physical laws such as balance and constitutive laws. 
As it will be showed herein below, balance law can be written in a discrete 
5 form by using the topological Equation (5) without approximation. The 
constitutive laws cannot be translated in algebraic form without approximation. 
Usually, they require the link between cochalns that belong to different cubical 
complexes. For example, in the case of dual complexes, two cochalns are 
linked by an algebraic linear system. This transformation is called "codual 

10 operation". More generally, 0-cochains represent local quantities. However, q- 
cochains (q>1) represent global quantities (e.g., an integral or the summation 
of a differential form) since they are associated with the algebraic structure of 
an edge, an area, a volume, etc. Thus, the cochain, coboundary, and codual 
are generic algebraic structures that can be instantiated by physical or 

15 mathematical laws. The exact translation of given problem in terms of q- 
cochalns is possible if one is able to find the basic laws that can be described 
without approximation by either cochalns or the topological Equation (5). 

In the previous example, the coboundary can be interpreted as follows: 
20 assuming that the vector field is conservative F = Av, ^AvjJs ==v(b)-v{a) 

constitutes a coboundary operator since it may be written as 
^0-^0(^1) = ^o(^i^i)» where a and b are the faces of Ci. Similarly, 
jdiv(y)dS^ jVMds, where n is the outward unit normal vector to Scj, 

constitutes a coboundary operator since it may be written as 
25 SiF^(c2) = i^i(92^2) • This example shows that the coboundary operator may be 
. an exact discrete representation of the fundamental calculus theorems (line 
integral and Gauss). 



Thanks to the chain, cochain, coboundary and codual concepts, the 
30 image model of the illustrative embodiments of the present invention can take 
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Into account the mathematical or physical laws related to the application. It can 
thus be Instantiated to the various problems of image processing and 
computer vision. The underiying computational framework Is strongly different 
form existing frameworks. For example, let us consider physics based 
problems such as optical flow, diffusion and deformation. Existing frameworks 
can be summarized as follow: 1) modeling by partial differential equation; 2) 
resolving the PDE by using numerical analysis scheme or Fourier space. 

The computational framework, according to the illustrative 
embodiments of the present Invention, can be summarized as follow: 1) 
identification of basic laws associated to the problem (Block 6701 of Figure 
6Z); 2) definition of an image support including the number of cubical 
complexes and the dimension of the cubes (e.g., for the case of a multi- 
resolution processing) (Block 6702 of Figure 67); 3) definition of global and 
local quantities (Block 6703 of Figure 67); 5) Instantiation of the coboundary 
and codual operators (Block 6704 of Figure 67); and 6) the resolution of 
resulting algebraic system (Block 6705 of Figure 67). The advantages of this 
framework are described hereinbelow and will be better defined in two 
practical examples. 

To summarize, the coboundary operator links quantities associated with 
the faces of an n-pixel. Codual operator links quantities associated to 
complexes of an image. If a given problem can be broken down into basic 
laws, the cochain and coboundary are the discrete representation of these 
basic laws. Cochain, coboundary, and codual are generic concepts that can be 
instantiated by physical or mathematical laws. Thus, they can be used in 
various computer vision and image processing problems. 

An illustrative example: 

Let us consider the linear isotropic diffusion In gray level images which 
is a physics-based problem. One can find all details In reference [4]. For the 
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sake of clarity and without loss of generality, the analysis will be limited to 
considering the 2D global differential equation for heat flow in a homogeneous 
medium. Let V be a 3D region with boundary S Inside the flow. The rate of 
heat flow across boundary S out of V Is given by: 

5 

/ / / cr(x,t)dy = " / / / v.(Avr(x,t))dy 

J J Jv J J Jv (6) 

where x is a spatial vector, t represents time, and A is a positive constant. 

10 To resolve this problem, the image support is defined by a continuous 

scalar field T, the temperature (i.e., gray level), two dual cubical complexes 
(i.e., two chains), and three cochains. If only one cubical complex is used, two 
different orientations are associated with each 1-face. To overcome this 
problem, two dual complexes (primary and secondary) are used (see Figure 

15 3). Concerning the use of three cochains, it has been pointed out in reference 
[4] that this EDP is formed by three basic physical laws. Each cochain is 
associated with one law. The first is the thermal tension law (also known as 
Pick's Law), which states that heat flows from regions of higher temperature 
towards regions of lower temperature. The direction of the gradient vr is the 

20 direction of the largest increase in temperature; the heat flows in the opposite 
direction. Formally, this law Is written as follows: 

flr(x,t)=-Vr(x,*) 

25 The primary cubical complex is the support for this balance law. The 

orientation plotted on this cubical complex is the direction of the path on which 
the integral is computed. Let us assume that the 0-cochain is the temperature 
T associated with 0-pixel Cq. A 1 -cochain G associated with 1 -pixel Ci is the 
global thermal tension transferred by the two 0-pixels that are the faces of CiJ 

30 Consequently, the topological equation is: 



wo 2004/015631 ^0 PCT/CA2003/001214 

26 



By the linearity of cochains, this topological equation is valid for all 1- 

5 chains. 

The second law, called the heat source law, concerns the net outflow 
of internal thernnai energy at the point x and time t. This is a balance law. It Is 
given by: 



10 



a(x,*) = V.g(x, t) 



(9) 



When V.q(x,t) > 0, the outflow is positive and thermal energy must 
flow away from x. Similarly, if V.q(x,t) < 0, the inflow is larger than the outflow 
15 and thermal energy increases at x. The equilibrium for a diffusion process is 
attained if V.q(x.t) = 0. 

Let us consider the secondary cubical complex. The orientation plotted 
on this cubical complex is the direction of the flow. The 2-cochain E 
20 associated with the 2-pixel cz is the global heat transferred by the faqes of Ca 



S(C2) = SiQ{C2) = Qid2C2) (^Qj 



This topological equation is valid for all 2-chains. The third law is 
25 constitutive (it depends on the medium feature). It makes the link between the 
flow density law and the heat source law. It is given by: 



q{x,t) =Xg{x,t) 



(11) 
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Thi6 equation cannot be discretized exactly, since it involves global 
quantities defined on two dual complexes. Consequently, the 2-cociiain ± 
cannot be computed without approximation from the 1-cochain G. For 
example, they can be linked as a linear equation system KQ = ± , where A is 
5 the coefficient matrix. Figure 4a gives an example of original image and Figure 
4b and example of image smoothed using this computational scheme. 

The data structure associated to the linear diffusion problem Is defined 
by: 1) two dual cubical complexes; 2) two cochains G and A for global 

10 quantities and a scalar field T; 3) two coboundary operations for balance laws 
and a codual operation tinat represents the constitutive law. The framework is 
summarized as follows: g is approximated by a bilinear polynomial. The 
cochain G is computed using a line integral, q is computed using the 
constitutive law in Equation (11). The cochain ± is computed using Gauss's 

15 theorem. Finally, a system of linear equations obtained from Equation (11) is 
obtained where the unknown variables are T. It should be noted that the 
cochains in Equations (8) and (10) are computed without approximation. 

The image model according to the illustrative embodiments of the 
20 present invention and described hereinabove may generally be characterized 
by three major points: 1) the image support and quantities are defined 
separately and then linked together via algebraic language; 2) the pixel is 
dimensional and is written in an algebraic form; 3) both local and global 
quantities are represented by the cochains and coboundary operators. 

25 

Each of these specificities will now be discussed to show their 
straightforward consequences for image processing. 

The separability of the image model allows a distinction between 
30 image variables and image quantities. The image variables offer numerous 
possibilities for existing mathematical formulations such as the use of 
algebraic topology to help in the design of algorithms. For example, binary 
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image algorithms are written as algebraic systems [16]. The well-defined 
quantities allow the use of physics, vector analysis or diffei^ential forms in the 
design of algorithms. Taking image support and image quantities together, 
well-known problems such as those of diffusion and optical flow in gray level 
5 images can be written as algebraic systems. Furthermore, the transfer of 
quantities between a given domain and its boundary is straightforward, using 
the concepts of cochain and coboundary as a general framework. For 
example, as shown hereinabove, in vector calculus, this transfer is easier 
thanks to the three fundamental calculus theorems, the line integral, Stokes's 
1 0 theorem, and Gauss's theorem. 

Unlike existing image models, by considering the pixel as dimensional 
primitive the connectivity paradox of the image support is avoided [8]. That is, 
the well-known Jordan theorem is fulfilled. Its decomposition into faces and the 

1 5 use of cubical structures such as chains make the dimension of the image 
explicit. Algorithms designed according to this formalism ^ operate in any 
dimension. This fact overcomes the traditional limitations that are faced in 
designing an algorithm, say in one dimension, extending it to two dimensions, 
and then to three dimensions, and so on. Each extension step may be a 

20 difficult task. 

The definition of the cochain depends on the problem that is dealt 
with. Thus, this image model offers real flexibility for the integration of 
mathematical objects (scalar, vector, tensor) and physical laws (balance, 

25 constitutive). Furthermore, the use of global quantities associated with an n- 
pixel implies noise reduction. In fact, global quantities are computed from the 
field by using the integral or the discrete summation. As the opposite operation 
from, the differentiation, which enhances high frequencies of the image, the 
integral performs a smoothing operation. It allows us to reduce the order of the 

30 derivative used in an image-processing scheme. Consequently, the 
introduction of global quantities may allow the use of higher-order derivatives. 
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Another contribution of the image model according to the illustrative 
embodiments of the present invention concerns the numerical scheme used to 
solve nonlinear problems such as the diffusion problem and elastic matching. 
It should be reminded that, usually, a problem is formulated and then a 
5 numerical analysis scheme is used. The numerical analysis scheme may not 
have been derived for exactly this formulation and many approximations must 
be made. The explanation of Intermediate results is not available. 
Consequently, no clear idea is available about the convergence of the 
numerical analysis scheme and the numerical results obtained may be a broad 

10 approximation of the desired solution. Based on the problems tackled, it is 
concluded that in the image model presented here, the numerical scheme Is 
deduced from the problem model with little or no approximation [4, 12]. In fact, 
various problems may be broken down into basic laws and then reformulated 
in terms of cochains and coboundaries. Thus they can be written as linear 

15 algebraic systems and solved. 

PRACTICAL EXAMPLE #1: PHYSICS-BASED RESOLUTION OF DIFFUSION 

AND OPTICAL FLOW 

An alternative to Partial Differential Equations (PDEs) will now be 
described for the solution of three problems in computer vision: linear isotropic 
diffusion, optical flow and nonlinear diffusion. These three problems are 
modeled using the heat transfer problem. Traditionally, the method for solving 
physics-based problems such as heat transfer is to discretize and solve a PDE 
by a purely mathematical process. Instead of using the PDE. the global heat 
problem can be decomposed into basic laws. It will be demonstrated that 
some of the basic laws admit an exact global version since they arise from 
conservation principles. It will also be showed that the assumptions made on 
the other basic laws can be made wisely, taking into consideration knowledge 
about the problem and the domain. The above-described image model will be 
used to allow encoding of physical laws by linking a global value on a domain 
with values on its boundary. The resulting algorithm performs In any 
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dimension, Tfie numerical sclieme is derived In a stralghtfonvard way from the 
problem modelled. It thus provides a physical explanation of e^ch solving step 
in the solution. 

5 Background 

in recent years, Partial Differential Equations (PDEs) have attracted 
increasing interest in the field of computer vision. Since PDEs have been the 
subject of much study by numericians, powerful numerical schemes have been 
10 developed to solve them. Consequently, domains such as image 
enhancement, restoration, multi-scale analysis and surface evolution all 
benefit greatly from PDEs [25]. 

One important class of equations governing certain physical processes 
15 . is the linear elliptic PDE of the general form known as the Helmholz equation: 

VM^) + Kx)zi(x) = /(x) ^^2) 

where x denotes a vector in the n-dimensional space, u(x) is the dependent 
20 variable, V^is the Laplacian operator, and p(x) and f(x) are spatial functions. 
When p(x) = 0. this corresponds to the Poisson equation [21 , 32] (also known 
as the non-homogeneous Laplace equation [30, 44]). One of the physical 
processes governed by Equation (12) is steady-state heat transfer. 

25 In the field of computer vision, Equation (12) may arise from two 

approaches. The first is variational calculus. As a matter of fact, many 
problems such as shape from shading [39], surface reconstmction [32, 40] and 
the computation of optical flow [29] can be formulated as variational problems. 
The solutions to these variational problems are given by Euler-Lagrange 

30 equations, which are in the form of Equation (12) [31, 32], The second 
approach is physics-based. For example, diffusion processes arise from heat 
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equations and shock filters from work in fluid mechanics [25]. For both the 
variational and the physics-based approaches, the resulting PDEs are 
continuous and have to be discretized. 

5 Traditionally, the discretization of PDEs in computer vision has been 

done by applying finite difference methods [23, 31, 39, 40]. Equation (12) is 
solved iteratively using either a direct Fourier-based Poisson solver for each 
iteration [39], finite elements [24], or spectral methods [32]. Iterative methods 
such as those in [39] do not ensure convergence unless smoothness is very 
10 high [21]. 

Existing methods for the resolution of problems involving PDEs can be 
summarized as follows: 1) identification of basic laws; 2) combination of the 
basic laws in oi-der to write the PDEs; 3) discretization of the PDEs; 4) 
1 5 resolution of the PDEs via a numerical method. 

This process, which has been used in various fields of application, is 
purely mathematical. Consequently, it has the following drawbacks: 1 ) Some 
quantities involved in the solution process do not have a physical 
20 interpretation; 2) This lack of interpretation is manifested in intermediate 
solutions involving iterative processes and since these solutions cannot be 
physically explained, discovery of the optimal solution cannot be ensured in an 
optimal time. 

25 Solution 

To overcome these drawbacks, an alternative to PDE resolution in the 
context of the heat transfer problem is proposed and will be described 
hereinbelow. 



Generally; basic laws in physics-based problems are combined into a 
global conservation equation [42] that is valid on the whole body or a part of it. 
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A local conservation equation (PDE) is then obtained by considering the 
particular case of a part of a body reduced to a point. 

In discrete problems such as those encountered in computer vision, 
5 the continuous domain is subdivided into many sub-domains in which there is 
only one value available, which can be considered as a global value. 
Therefore, instead of using the PDE, this illustrative embodiment of the 
present invention proposes to use the global conservation equation directly on 
each sub-domain. 

10 

In order to handle these physical laws which link global values at 
points, lines, surfaces, volumes, etc. the image model with roots in 
computational algebraic topology described hereinabove is used. This model 
makes it possible to represent global values on entities of any dimension at the 
15 same time. 

The above described methodology presents a number of advantages: 

1) Many of the basic laws arise out of conservation principles and hence they 
20 are valid either at a point (local form as in Equation (12)) or over an entire 
region (global form). Fundamental theorems of calculus such as the 
Gauss, Green and Line Integral theorems allow the computation of the 
coboundary operator without any approximation. 

25 2) Some laws require approximations that can be performed wisely, taking 
into account knowledge about the problem and the domain. 

3) The intermediate results have a physical explanation because they 
represent physical quantities. For that reason, every step has a physical 
30 interpretation. Thus there are no longer problems of non-optimality of the 
solution, because we avoid non-temporal iterative methods. 
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4) As mentioned earlier, this metliod can be used with other physics based 
problems by applying the appropriate basic laws [36]. 

5) Thanks to the image model, the resulting algorithm performs in any 
5 dimension. 

6) The computational rules associated with the coboundary operator can be 
changed without changing the formalism of the operation itself. 

10 7) The same formalism can be used for pixel-based and other types of 
decomposition of the Image (e.g. regions). 

In order to validate the method, the equation for steady state heat 
transfer is resolved in two applications: the linear diffusion and optical flow. 
15 These problems generate equations of the form of Equation (12) or its global 
version. The present methodology can also be used to resolve unsteady heat 
transfer with no source and we apply it to non-linear diffusion. 

Physical Principles (explanation of the physical foundations of the heat transfer 
20 problem) 

Two interesting particular cases for diffusion and optical flow problems 
can be distinguished: steady-state heat transfer and unsteady heat transfer 
with no source. 

25 

Two important classes of laws are present: conservation and 
constitutive laws. Conservation laws are independent of the properties of the 
material, whereas constitutive laws are specific to them. The physical 
properties associated with a moving object are energy, work and heat. In what 
30 follows, each of these quantities are described. 



Energy Modeling 
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Some quantities for a continuous body occupying a volunne V 
bordered by a surface S in a 3D space will first be defined. Such a body can 
be said as composed of an infinite number of particles (as many as desired), 
these particles being the smallest elements [33]. Figure 5a illustrates such a 
body. At time t, a particle labelled X occupies a specific position: 



Each particle can move in space, so a velocity vector is associated 
with it at time t: 

v*(X,t) = ^=v(x,*) 

Physical quantities can be associated with a particle labelled X 
(material description) or a position x in space (spatial description). For 
example, v*(X,t) is the material description of the velocity of particle X and 
v(x,t) is the spatial description of the particle located at position x. For the 
present purpose, spatial descriptions are used to derive the heat transfer 
equation. Vector n(x,t) is the outward direction of the surface at point x. 

The mass Am of a small amount of volume AV of a body is a 
measure of its inertia (tendency to resist motion). Tlie term mass density, p is 
used to denote the following quantity: 



x = (a:(t),y(t),-2?W) 




Am 

nm — 

AV-*0 AV 



p (x,t) is thus the mass density of the particle located at x at time t. 
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Two kinds of energy are associated with a moving object: kinetic and 
internal energy. Kinetic energy Is a measure of the state of motion of a body: 
the faster the body moves, the greater Its kinetic energy [28]. Because It is a 
measure of inertia, kinetic energy also takes the mass into account. For a 
5 particle located at x at time t, the kinetic energy Is thus defined as 



where is the dot product. To obtain the kinetic energy for the entire body at 
10 time t, K(x,t) Is integrated over the volume V: 



The hotter the body, the greater its internal energy. At time t, each particle has 
an internal energy density s{x,t) associated with it. The intemal energy 
density is proportionai to the temperature of the particle T(x,t) with a material- 
specific heat constant c; that is, e (x,t) = cT(x,t). For the entire body, the total 
20 internal energy is integrated over the volume V: 



= ip(x, t) (v(x, t) • v(x, t)) 



^ii) = 1/ ff^Pi^^i) (v(x,t) • v(x,«)) dV 



where dV is an infinitesimal amount of the volume V. 



Intemal energy is a measure of the state of temperature of a body. 




(13) 



The total energy for the body can now be defined as K(t) + E(t). 



Work Modeling 
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Let us suppose that a body is submitted to an external force fe(x,t) 
(e.g. a traction force) and an internal density force b(x,t) (e.g. gravity). Figure 
5b presents the action of external and Internal forces. Work is defined as the 
energy transferred to a body by means of a force acting on the body [28]. 
5 Work is negative when the energy is transferred from the body. Suppose that 
F(x) is an Internal or external force that is constant over tinne, acting on a 
particle located at x during an amount of time t. This force will produce a 
displacement of the particle to position xi. This displacement Is Ax = xi - x. 
The work w(F, x) done by this force during this time is: 

10 

i£;(F,x) = F(x)- Ax 
and the instantaneous power P(x, t) is: 

F(x,i) = lim ^^^^ = F(x) . liin ^ = F(x) • $ = F(x) - v(x,t) 

External forces act essentially on the surface of the body. The 
instantaneous work Pe(t) done by the extemal forces on the entire body is thus 
the result of the integration of the extemal power over the surface: 

20 

where dS is an infinitesimal part of the surface of the body. 

25 Defining b(x, t) as the internal density force, the internal force is thus 

p (x, t)b(x, t). The rate of work over time done by the internal forces on the 
entire body is obtained by integrating intemal work over the volume: 
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Pi(t) = fjf^ *) (b(x, t) • V(X, t)) dV 

The total work is thus P(t) = Pe(t) + P|(t) 

5 Heat Modeling 

Heat can be defined as the energy transferred to a body owing to a 
difference in temperature- The heat flow density vector q(x, t) is a measure of 
the rate of heat conducted into the body per unit area per unit of time. How 
10 q(x, t) is defined will be explained later. The external heat addition rate over 
time is the amount of heat coming from outside the body and entering by its 
surface.. It is computed by projecting q(x, t) onto the inward normal vector 
(-n(x, t)) and integrating this projection over the surface: 



Qe{t) = ff q(x, t) • (-n(x, t)) dS 
15 ''''S (14) 



Now if a body has a rate of heat generation per unit of volume, and 
time r(x, t), the internal rate of heat addition over time is computed by 
integrating r(x, t) over the volume: 



20 



Qi{t) = j ff^Pi^**) ^(X'*) ^ 



Figure 5c shows q(x, t) and r(x, t). To simplify, the source o- (x, t) is 
defined as the rate of heat generated in a particle located at x per unit of 
25 volume and time: 
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In many cases, this source is known. However, it couid also be a 
linear function of the temperature: a (x, t) = a(x, t) + b(x, t)T(x, t) [38]. The 
total rate of heat addition over time is thus: 

5 

Q(*) = Qe(i)+Qi(*) 

Energy Conservation Law 

.10 A class of equations in continuum mechanics are those describing the 

conservation (equilibrium) principles. They express the conservation of certain 
pliysical quantities (mass, momentum, energy, etc.) over an entire body and 
as such they take the form of global equations over the whole body or a part of 
it [33]. , 

15 

Conservation principles can be seen intuitively as follows: the change 
in the total amount of a physical quantity Inside a body is equal to the amount 
of this quantity entering or leaving the body (through the boundary) and the 
amount generated or absorbed within the body. These laws are applicable for 
20 all continuous, materials, moving and stationary, deformable and non- 
deformable, and must always be satisfied. The global conservation equations 
can then be used to derive their local counterparts, called the associated field 
equations, which are valid at each point of the body including its borders. 

25 The first law of thermodynamics, vvhich Is relevant for the 

understanding of the heat transfer equation will now be discussed. This law 
involves both kinetic and internal energies and states that the total variation of 
energy in a body (or a part of a body) is the result of the time rate of work and 
the rate of heat addition combined: 



30 
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^ (16) 

For heat transfer, the only interest resides in the case of immobile 
bodies, that Is v = (0,0.0). n(x.t) = n(x) and p{x.t) = pix). Equation (16) thus 
5 becomes: 



which now states that the thermal energy variation in a body Is due to internal 
10 heat production added to the heat flowing Into the body, U^ing the divergence 
tlieorem for Qe [33] and recalling that ^ (x,t) = cT(x,t), we obtain the themial 
energy conservation law: 

15 (17) 
where V. is the divergence operator. To simplify, let us define the temperature 

variation h(x,t) = -^T(x.t). Equation (17) is a conservation equation and Is 

at 

thus valid over the entire body, a part or a point of this body. Consequently, 
the integral signs can be taken off: 

20 

cp(x) fe(x,^) = -V •q(x,t) + <7(x,<). 

Tbemudeii«:^vailatioii Rate of heat eatenng Rate of heat genendioa 

This equation Is said to be local, whereas Equation (17)' is said to be 
global. The themnal energy variation Is called the unsteady term, the rate of 
25 heat entering Is called the diffusion term and the rate of heat generation Is 
called the source term. 



15 
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Based on Equation (18), two cases are be considered: tlie steady 
state case and the unsteady case witli no source. The term steady simply 
means that there Is no variation of the thermal energy of the system over time. 
5 That is, the left side of Equation (1 8) is null: 



This implies that the heat diffusion compensates for the internal heat 
10 production: 



In the unsteady case with no source we have: 
cr(x,t) =0 



(21) 



which means that the time variation of thermal energy is explained by the heat 
diffusion alone: 



cp(x)hiK, t) = - V • q(x, t) ^^^^ 



Constitutive Principles 

25 In Equation (18), there are three unknown variables: /7(x). h(x,t) and 

q(x,t). Let's look at the example of q(x,t). Suppose that we can measure the 
time variation of the thermal energy (left side of the equation) and also of the 
temperature T. We know that q(x,t) is related to the temperature, but since 
different materials usually have different diffusion properties, the missing 
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equation q(x,t) = f(T,x,t) must depend on properties of the material we are 
studying, such as its homogeneity and type of diffusivity. Consequently, the 
system of equations contains more unknown variables than equations and the 
function f(T,x,t) must be added to the system formed by Equation (18). This is 
5 due to the difference in material properties. Different materials behave 
differently, but are subject to the same conservation laws. Constitutive 
equations such as f(T,x,t), which reflect the internal constitution of materials, 
allow us to complete the system of equations. 

1 0 Decomposition into Basic Laws 

As Indicated in the foregoing description, conservation equations are 
always valid regardless of the materials, whereas constitutive equationsi are 
dependent on their properties. When solving directly PDEs like Equation (18) 

15 in a discretized context with methods such as the finite differences approach, 
one makes global assumptions about the time and space behavior of the 
diffusion, energy variation and source terms without taking into account the 
nature of the basic laws underlying the problem. Some of these do not require 
any approximation since they come from conservation principles. Also, a more 

20 physically realistic solution can be obtained by choosing a proper 
approximation for each basic law arising from a constitutive principle. 
Consequently, we propose to decompose the terms of Equation (18) into basic 
laws. This equation can be brol<en down into one constitutive and two 
conservative laws for the steady state case. In the unsteady case, an 

25 additional constitutive and another conservative law must also be considered. 
Note that since the source term is often known, it will not be decomposed. 
Recalling that the diffusion term a (x,t) is the rate of heat entering the particle 
located at x at time t, then: 



30 



a(x,t) =-V -q(x,t) 



(23) 
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is a first basic conservation law. 

The second conservation law concerns the thermal tension. We first 
define the thermal tension vector g(x,t) as the vector representing the direction 
5 and magnitude of the greatest temperature decrease at a fixed time t. As g(x,t) 
is source-oriented (from hot to cold), a minus sign must be Inserted before 
VT(x,t) which represents the direction and magnitude of the greatest 
temperature increase: 

10 o\ w \ > / ^24) 

^ This equation is a second basic law. Since the thermal tension is the gradient 
of a scalar field, it is by definition a conservative field in space, it can be said 
that -T(x,t) is the potential field of g(x.t) [26, 41]. 

•15 

The third law is a constitutive law. The heat flow density q{x,t) is 
defined as the quantity and the direction of the heat flowing into the particle 
located at point x at time t. It is represented by a vector and greatly depends 
on the behavior of the material, in the case of uniform, homogeneous 
20 materials, it has been proven experimentally by Fourier [20, 27] that q(x,t) is 
directly proportional to the difference of temperature relative to neighbors of 
this particle: 



25 



q(x,i) =Ag(x,*) ^25) 



where ^ is a material-specific thermal conductivity constant. The value of A. is 
known for many types of materials. Equation (25) is called the Fourier heat 
conduction law. For a non-homogeneous material, we consider that it has the 
behavior of a homogeneous material on an infinitesimal patch, but the 
30 conductivity changes with each patch; that is: 



wo 2004/015631 




PCT/CA2003/001214 



q(x,*) = A(x,i)g(x,<) 
For the unsteady case, the fourth basic law is: 



(26) 



where s (x,t) is the thermal energy variation (the unsteady term). This equation 
is a constitutive one because it involves p(x), whicfi is material-dependent. 

10 Finally, the fifth basic law is related to the temperature variation and is 

a conservation law: 

Mx,*) = ^r(x,i) 

at (27) 

Considering only the temperature Tx(t) of the particle located at x, we 
15 reduce the basic law to a 1 -dimensional equation and we can thus say that 
Tx(t) is a conservative field in time-space. 

To summarize, three basic laws for the diffusion term of equation (18) 
have been defined, that is: 

20 

a(x,t) = -V.q(x:,t) 
q(x,e) = Ag(x,i) 
g(x,i) = ~VT(x,t) 

25 There are also two additional basic laws for the unsteady temi, that is: 



e(x, t) = cp(x) /i(x, t) 
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Combining all these elements, the following relation is obtained: 
cp(x,*)^T(x,t) = V . (AVT(x,i)) + <r(x,<) 

(28) 

Discrete Representation of Images 

Some algebraic tools used to model images will now be recalled from 
the above description. An image Is composed of two distinctive parts: the 
image support (pixels) and some field quantity associated with each pixel. 
This quantity may be scalar (e.g. gray level), vectorial (e.g. color, multispectral, 
optical flow) or tensorial (e.g. Hessian). The image support is modelled In 
terms of cubical complexes, chains and boundaries as described in the 
foregoing description. With these concepts, it Is possible to give a formal 
description of an image support of any dimension. For quantities, the concept 
of cochains has been Introduced, these cochains being representations of 
fields over a cubical complex. For the use of these concepts in Image 
processing, see [16]. 

As discussed hereinabove, an image is a complex of unit cubes 
usually called pixels. A pixel / c SR" is a product : 

7 = 7i X /a X . . . X /„ 

where Ij is either a singleton or an Interval of unit length with Integer endpoints. 
Thus Ij is either the singleton {k} and is said to be a degenerate inten/al. or the 
closed interval [k. k+1] for some A:eZ. The number q e{0.1 n} of non- 
degenerate Intervals is by definition the dimension of which is called a q- 
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pixel. Figures 6a-6c iilustrate tiiree elementary pixels in SR^ . For ^ > 1 , let 
*/ = {fco5^i>- -^^-i} ordered subset {1,2. ...,n} of indices for which 

Ij^ = [ajybj] is non-degenerate. Let us define: 

and 

Bk^a = /ix.../fc^.^ix{6j}x Jfc^+ix ...x/n 

The Af^^ and the -B^fe^are called the (q-l)-faces of cr. One can define 

10 the (q-2)-faces, .... down to the 0-faces of a in the same way. The faces of / 
different from y itself are called its proper faces. 

By definition, a natural orientation of the cube is assumed for each 
pixel. Suppose that / denotes a particular positively oriented q-pixel. It is 
15 natural to denote the same pixel with opposite orientation by -y. Examples of 

orientations are given in Figures Ba-Sb. A cubical complex in 91" is a finite 
collection K of q-pixels such that every face of any pixel of the image support 
is also a pixel in K and the intersection of any two pixels of K is either empty or 
a face of each of them. For example, traditional 2D image models only 
20 considered pixels as 2D square elements. The definitions presented above 
allow us to consider 2-pixels (square elements), 1 -pixels (line elements) and 0- 
pixels (point elements) simultaneously. 

In order to write the image support in algebraic form, the concept of 
25 chains is introduced. Any set of oriented q-p!xels of a cubical complex can be 
written in algebraic form by attributing* to them the coefficient 0,1 or -1, if they 
are not in the set or if they should or should not be taken with positive 
orientation, respectively. In order to represent weighted domains, arbitrary 
integer multiplicity is allowed for each q-pixel. 
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Given a topological space Xam" in ternns of a cubical complex, we 
get a free abelian group Cq(X) generated by all the q-pixels. The elements of 
this group are called q-chains and they are formal linear combinations of q- 
6 pixels [16]. A formal expression for a q-chain Is Cg^^^^^/i where 

The last step needed for the description of the image plane is the 
introduction of the concept of a boundary of a chain. Given a q-pixel / , we 

10 define its boundary as the (q-1)-chain corresponding to the alternating sum 
of Its (q-l)-faces. The sum is taken according to the orientation of the (q-1)- 
faces with respect to the orientation of the q-pixel. A (q-l)-face of y is said to 
be positively oriented relative to the orientation of / If its orientation is 
compatible with the orientation of / . By linearity, the extension of the definition 

15 of boundary to arbitrary q-chains Is easy. For instance, in Figures 6b and 6c, 
the boundary of the 1 -pixel a is Xa - Xi and the boundary of the 2-pixel A is 
a + b-c-d; then a and b are said to be positively oriented with respect to the 
orientation of A but c and d are said to be negatively oriented with respect to 
the orientation of A. Let us notice that the boundary of a 1 -pixel is always the 

20 difference between Its boundary points. The boundary can be defined 
recursively. Given a (q-1)-chain and a q-chaIn defined as/g =rg-i^[a>b], 
the boundary of can be recursively written as: 

dj, = ^7,-1 X [a, b] + (-l)(«-^)(7^^, X {6} - 7^,, X {a})^^^^ 

25 

In order to model the pixel quantity over the image plane, an 
application F that associates a global quantity with all q-pixels y of a cubical 
complex is determined and is denoted by <F,r>. This quantity may be any 
mathematical entity such as a scalar, a vector, etc. For two adjacent q-pixels 
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n and F rn^st satisfy <F,Ji^ri'^^r2 >=^K <F>ri>+^ <F,r2 >» which 
means that the sum of the quantity over each pixel is equai to the quantity over 
the two pixels. The resulting transformation F:C^(X)-^9i is called a q- 

cochain and is used as a representation of a quantity over the cubical 
5 complex. 

Finally, an operator is needed to associate a global quantity with the 
(q+1)-pixels according to the global quantities given on their q-faces. Given a 
q-cochain F, we define an operator d , called the coboundary operator, which 
10 transforms F into a (q+1)-cochain SFsuch that: 

^ ' (30) 

for all (q+1)-chains The coboundary is defined as the signed sum of the 
15 physical quantities associated with the q-faces of /. The sum is taken 
according to the relative orientation of the q-faces of the (q+1)-pixeis of y with 
respect to the orientation of the pixels. Figure 7 presents an example of the 
coboundary operator for a 2-pixel. 

20 With this image model in hand, the basic laws described hereinabove 

will be used to rewrite the global heat transfer equation in algebraic terms [43]. 

Representation of the Heat Transfer Equation 

25 The process for representing the heat transfer equation in terms of 

algebraic topology can be summarized as follows. The Image support is 
subdivided into cubical complexes. Basic laws are applied to pixels of various 
dimensions. These laws involve the computation of global quantities on pixels, 
expressed as cochains. Some of these laws link global quantities on pixels 

30 with global quantities on their boundaries and hence are expressed as 
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coboundaries. The other laws are expressed as linear transformations 
between pairs of cochains. The topological fonnalism of cochain and 
coboundary is a generic one; that is, it does not offer computational rules. The. 
cochains must be instantiated depending on the problem to be considered. 

5 

The basic laws presented hereinabove will be reformulated in a 
topological way and then to give computational rules for cochains in the 
context of the heat transfer problem. Since we want to represent two kinds of 
global values over the spatio-temporal image, two complexes will be used. The 
10 first complex Is associated with global values corresponding to projections 
onto the tangential part of the domain (e.g. global thermal tension) while the 
second complex refers to values related to projections onto the normal part of 
the domain (e.g. heat entering a particle). These two distinct orientations (see 
Figures 8a and 8b give rise to two different complexes. 

15 

Global Heat Transfer 

Let us assume that an image has n spatial dimensions and r pixels. 
Suppose also that a time interval [to, tj can be split into i equal sub-intervals [tk, 
20 tk+i], k e [0, i-1]. Let us consider an n-complex representing the subdivided 

spatial support of the image K' . One can consider an (n+1).complex 
representing the spatio-temporal support of the image: 

= K'' X ltk,iM],Vk £[0,1- 1] 

25 

Now. let us consider r^, an (n+1)-pixel of 1^, and the following 
cochain, defined as <s,yg>. Thus, we therefore need to define which value 
to use as cochain s in the heat transfer problem. Let us define e as the 
global energy variation of the (n+1)-pixel : 
30 • 
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< IB >= / C(x, t) djB 

J^E (31) 



where dy^ is an Infinitesimal part of the domain represented by Ye ■ Now, 
using the global version of Equation (18), we obtain : 

5 

/ e(x, *) d^E = I *) djB + / ^(x, t) d^y^ 

J^B J^B J'TB 



From this equation, we define two more cochains, representing first, the global 
diffusion: 

10 



Yb 

15 



'B 

rfB (32) 

and second, the global source: 

iS,7jb >= / a(x,t)d7j5 
Thus, the following relation is obtained between the three cochains: 
<S,^B >=< V,^E > + <S,^E > 

(33) 

20 

The rules used for cochains ^and D are then decomposed into basic 
laws. The rule for cochain S is not decomposed since it is assumed that its 
global value is known on y^. Let us finally mention that both steady and 
unsteady heat transfer problems can be considered using Equation (33) by 
25 setting respectively, cochains € and S, to zero. 



1 
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Global Temperature Variation 

Let us consider another n-complex, , representing the subdivided 
spatial domain of the Image. An (n+1)-complex representing the spatio- 
temporal image can then be defined as: 

KF = KF'>< [tk, tk+i],Vk £[0,1- 1] 
Let us consider rir, a 1 -pixel of iST" defined as xix[tk,tk+i], ie[1,r|, 

f 

l<e[0,l-1] wliere Xi is a 0-pixel of K" . Let us also consider a 0-cochain T and a 
1-cocliain H sucli that: 

< W, 7H >=< ST, 7F >=< r, djH > 

(34) 

Figures 9a and 9b present examples of cochains T and H for of 
dimension 3. ' 

Applying Equation (29), It is found that the boundary of rn >s 5y^ = 
xix{tk+i} - Xix{tk}. According to the linearity of the cochain, the computational 
rule relating the global value associated to with the values at its boundary 
Xix{tk} and xix{tk+i} is: 

< Tydjff >=< T,JCi X {tfe+i} -Xi X {tf,} >=< r,Xi X {4+1} > - < r,Xi X {tk} 

(35) 

This equation is general and applies to many problems. To define 
which values to use as 0-cochain and 1 -cochain, let us take the global version 
of Equation (27) on and apply the fundamental calculus theorem: 
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Looking at this equation, It can be seen tliat it is similar to Equation 
(35). Tlius we define T = T(x. t). The location of the unl<nown temperatures to 
compute will correspond to the 0-pixels of . In order to fulfill Equation (34), 
the following relation is used: 



this relation being called the global temperature variation. These three 
equations are extended by linearity to a 1 -chain of iiT^ defined as ;^x[tk. tk+i], 



where ;^ls an arbitrary 0-chain of . 

Global Energy Variation 

We want to link cochains H and representing the global 
temperature variation and the global energy variation, respectively. For this 
purpose, a representation of Equation (26) is needed. The two cochains are 
not from the same cubical complex (H is from and s Is from K"" ), and 
moreover, Equation (26) is material-dependent; therefore they cannot be 
linked exactly. However, we can express this link as a linear transfomnatlon: 




(36) 



H 



Recalling Equation (31), the following relation Is obtained: 



< S,^B >= / c(x,t)d7B = / cp(x)/i(x,t) djB 





(37) 
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Unfortunately, the value of p{x) or h(x, t) is not known at all points of 
the volume. Consequently, these two fields are approximated over the volume. 
The approximations are denoted by p(x) and ^(x, t). For one 1 -pixel y^, 
5 defined as Xix[tk,tk+i], the approximation is performed piecewlse such that 
A (x. t) must fulfill Equation (36). 



(38) 



1 0 Equation (37) thus becomes: 



< >= / cp(x)^(x, t) djE = /e(c, n)=T 



(39) 



where dV is an infinitesimal part of which depends on the choice of the 
15 approximation functions p{x) and h(x, t) and the position of K"" with respect 

to i:^ . 

Global Diffusion 

20 Let us consider an n-cochain Q and an (n+1)-cochain D defined by the 

coboundary: 



< V,^B >=< 5Q,^E >=< Q.djE > 



25 Figure 10 presents examples of Q and D for iST^of dimension 3. Let us 

assume that the n-faces y^ of y^ are positively oriented relative to y^. 
According to the linearity of the cochain, the computational rule is: 
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<2?,7£>= ^ <Q,7Q,> 



Again, this equation is general; hence a global value is found for the 
5 (n+1)-cochain D, which can be computed by summing the global values at the 
boundary of rs • According to Equation (32), the following relation is obtained: 



10 The divergence theorem is applied to this equation to obtain: 



1 Js^\ 



where n(x, t) is the normal vector to an infinitesimal part of the domain 
15 represented by yg . This last equation is in the form of a coboundary (Equation 
(41 )), from the following relation is defined: 

< Q> iQi. >= / -q(x, t) • n(x, t) djQ, 

^rc^i (42) 

20 Again, the previous definitions can be extended by linearity to arbitrary 

(n+1)-chains of K' . And there Is absolutely no approximation in these 
equations. 



Global Thermal Tension 

25 

Let us consider a 1 -pixel of defined as yxtk, *€[0,/-l]. 
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where y is a 1 -pixel of whose boundary Is defined as 5y= Xj - xi , 
/,7e[l,r]. Let us also consider a 0-cochain T and a 1-cochain G defined by 
the coboundary: 

<g,jG >=< (JT, 7G >=< r, d-ya > ^^^^ 

Figures 11a and 11b present examples of cochains T and G for one 3- 
pixel of K" . The boundary of ra 's d/c = Xj x{tk} - x, x{tk}. According to the 
linearity of the cochain, the computational rule relating the global value 
associated with j^^to the values at xi x{tit} and 3q x{tij Is: 

< T, djG >=< T, X {tjfe} - Xi X {<!;} >=< r,xj X {tk} > - < r,Xi X {tk} > 

(44) 

To define which values to use as cochains G and T let us take the 
global fomn of Equation (24) on = 



/ g(x,t) 'dhfG = f ^ e{x,tk) • d7 = f ^ -VT(x,4) • <h 



(45) 



where dy is an infinitesimal part of y . Since g(x, t) is a spatial conservative 
field, we can apply the line Integral theorem [26, 41] saying that for a 
conservative field F(x) = Vf(x) and two points A and B, in an open connected 
region containing F(x). the Integral of the tangential part of F(x) along the 
curve R joining A and B is independent of the path (Figure 12): 

r F(x) dRi= r F(x) . di22 = r F(x) . dR^ = f(B) - /(A) 

J A Ja J a 
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From this theorem, Equation (45) can be rewritten as: 

(46) 

5 

Looking at Equation (46), it can be seen that it is similar to Equation 
(44). Thus T = T(x, t) Is defined. Consequently, the location of the unknown 
temperatures to be computed correspond to the 0-pixels of K'' which Is 
coherent with the conclusions hereinabove. In order to fulfill Equation (43), we 
10 have: 

< 6, IG >= / -g(x,i) • djG 

J-to (47) 

The previous definitions are extended by linearity to 1 -chains of K' 

t 

15 defined as r x{tk}, where / is an arbitrary 1 -chain of . 
Heat Flow Density 

The coboundaries <Q,dyE> (Equation (40)) and <r,dr^> 
20 (Equation (43)) provide exact global versions of Equation (23) on and 
Equation (24) on K'' , respectively. In order to complete the diffusion term. 
Equation 25, which links local values g(x,t) and q(x,t) Is represented. Equation 
(25) is a constitutive equation and cannot be represented by a topological 
equation. However, a relation transforming cochain G into cochain Q can be 
25 found: 



<Q.Jg> ^<Q,7q> 
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as a global counterpart for Equation (25). To find this transformation. Equation 
42 is recalled: 

< Q, jQi >= / -q(x, t) • n(x, t) dqq, = / -Ag(3c, i) • n(x, *) dr^Q^ 

5 

this equation relating cochain Q to field g(x,t). Unfortunately, field g(x,t) Is not 
known, so that this equation has to be approximated with a field g (x,t). Let us 

consider r„ . an n-pixel of K" defined as r„=rx^ i^} . A: e [O,/ -1] where is 

an n-plxel of K'' . This approximation Is performed plecewise such that for 
10 each 1-face y^of y„ , g (x,t) satisfies: 



/ -g(x, • dij =< TG > 



(48) 



where dR Is an infinitesimal part of the domain represented by /f,. Equation 
15 (25) is then applied to obtain q (x,t): 

q(x,*) = Ag(x,*) 
at all points of the domain. Equation (42) becomes: 



20 



< Q>7Qi >= / -q(x,*) • n(x,i) djQ, = fg(\ G) 

(49) 



The transformation that is looked for is thus: 



25 



A = /^(A,g?) 
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which depends on the choice of an approximation function g (x,t) and the 
position of AT* with respect to . 

5 Boundary Conditions 

The decomposition process that has been presented herein is carried 
out with the assumption that all the needed quantities surrounding a pixel are 
Icnown. For Instance, in the steady state heat transfer problem, for a particular 
10 (n+1)-pixel, the cochain S Is known for all other surrounding (n+1 )-pixels, that 
is, there are as many equations as variables. 

Unfortunately, this assumption is not verified at the borders of the 
image. Thus, as in solving the PDE, certain boundary conditions are imposed 
15 to specify the gray-level conditions at the boundary of the image. For instance, 
these conditions may prescribe the values of either cochain T (Dirichlet 
boundary conditions) or cochain Q (Neumann boundary conditions). 

. Summary of the Algorithm 

20 

The algorithm used to find an expression of the temperatures at time 
tk+1 as a function of the temperatures at time tk will now be summarized. The 
Input data for this algorithm are the cochain S and the Dirichlet boundary 
conditions. That is, T is known for all pixels on the boundary of the image, 
25 which includes the values at time to. 

1 . Choose the positions for and . 

2. Compute ^ as a function of H: 

(a) Choose the approximation functions ^(x.t) and p{x). 

(b) Apply Equation (38) to find h (x.tk, t,) as a function of H. 

30 (c) Apply Equation (39) to find the transfonmation r , expressing ^ as a 
function of H. 
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10 



20 



25 



3. Apply r and Equation (34) to find ^ as a function of T. 

4. Compute Q as a function of G: 

(a) Choose the approximation function g (x,t). 

(b) Apply Equation (48) to find g (x,t) as a function of G. 

(c) Apply Equation (49) to find the transfomiatlon A, expressing Q as a 
function of G. 

5. Apply Equation (40), A and Equation (43) to find D as a function of T. 

6. Apply Equation (33) to obtain an equation of the temperatures at time tk+i 
as a function of the temperatures at time tk 



Figure 13 presents an overview of the computational scheme. 

Applications 

1 5 Linear Isotropic Diffusion 

One of the most direct applications of the heat transfer equation is the 
isotropic diffusion of gray-level intensities; that is, smoothing. For a 2D image 
l(x), with X = (x,y), the resolution of the PDE: 



^J(x,t)==V^/(x,t) 



is equivalent to the convolution: 

I{x,t) = {I*gt)(pc) 

where 



(50) 



1 _/'«!±2i> 
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is a Gaussian with variance cr^ =2^ [25]. One can see t as tlie scale of tiie 
smootliing operation. Let us assume that the Laplacian image at scale t: 



5 ii(x,t)=V^/(x,t) ^5^^ 

is known. One can consider this equation as a steady state heat transfer 
problem with T(x,t) = l(x,t). a (x,t) = ~L(x,t) and A = 1 . 

10 It is desired to solve Equation (51) for local l(x,t) located at the center 

of each image pixel. Employing the process presented hereinabove, we first 
position the two cubical complexes representing two subdivisions of the image 

plane. As stated hereinabove, the primary complex is defined with 0-pixels 

corresponding to pixel centers. For the sake of simplicity, K' corresponds to 
15 the image pixels; that is, the secondary 2-pixels are rectangular and 

symmetrically staggered relative to the 1 -pixels of and the 1 -pixels of 

K'' intersect orthogonally in the centers of the primary 1 -pixels. Since there is 
no variation in steady-state heat transfer over time, the time parameter is 

dropped. This means that ^K^ , K' =K'' and the time integral in cochain 
20 computation are dropped. It can be seen that the approximation function fg 
depends on the position of iiT^with respect to . 

Figure 14 shows the two complexes for a 5 x 5 image. Positioning the 
1 -faces of such as each passes through the center point between two 0- 
25 pixels of allows us to compute a polynomial function of order 1 with the 
same accuracy as that obtained using one of order 2 [37, 34], 



A global value for the 2-cochain <S,y^ > is needed. If it is assumed 
that a pixel value represents the global value of intensity, < S^y^ > = -L(x) can 
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be directly set. This assumption is reasonable if image acquisition is 
considered as a process which accumulates the total number of photons within 
a global area corresponding to the pixel [22]. 

An approximation function g (x) is chosen. For simplidty, we assume 
that g (x) arises from a bilinear approximation, that is: 

i(x) = (a + 62/) • z + (c+ die) • J 

Given a 2-pixel y^of K" , g{x) satisfies Equation 48 for each 1-face 
of . As an example, let us find the coefficients a, b, c and d for such a pixel 
defined as in Figure 15: 



ft = y 


— g(a:, 0) • tdx 


^ J 


( -g(A,2/)-/(% 
0 


= \ 


^ -g(a;, A) • idx 
0 




' -S{0,y)'jdy 



0 



from which 

(52) 

is obtained, g (x) Is thus a piecewise function of G. but as G is computed from 
T, and g (x) can also be expressed as a function of T. For each primary 2- 
pixel, Equation 25 is applied to obtain q(x)= g (x). 
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The next step is to compute <Q,rQ >^or K' from Equation (49). Each 
secondary 2-plxel intersects with four primary 2-pixels, y^^ , y^^ , y^^ and 
y^a . There are four segments in the approximation function q (x) 
corresponding to the four primary 2-pixels; that is, §'„(x), §'^(x), §;(x) and 
?rf(x). Figures 16a-16c illustrate y^. Cochain <Q,yQ > corresponding to the 
four 1 -faces of y^ Is found by: 



Qi = -qa( A/2, y) 'tdy+ -qft( A/2, y) • idy 

Jo J'-AJ^ 



10 4 8 8 



^2 = / -q6(a?,-A/2)-(-^i)da7+ / ^q^{x,-A/2) ^ {^j)dx 

Jo J^A/2 

4 8 8 

Q3 = I -qa(a;,A/2) / -qrf(x, A/2) • jrfx 

^^O J-A/2 

^ 3g4 g2 gll 

4 8 8 



/•A/2 /.O 

^ = / -qrf(- A/2, 2/) • + / -qc(- A/2, y) • (^Qa?/ 

^'O J'^A/2 



15 ^~A/2 
3^10 ^ ^12 _ 

~ * ^ « (53, 
Using Equation 41 , we obtain: 



20 



< 2>,7E >= Q1+Q2 + Q3 + Q4 



(54) 
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Substituting Equation (43) in Equation (53), Equation (53) in Equation 
(54) and Equation (54) in Equation (33), kS^Ye > can now be expressed as a 
function of T. As an example, <S,rs> is presented for a 2-plxel y^, and 
5 defined as in Figure 16a-16c: 

< 5, 7* >= -375,0 + \ [75.1 + 71.0 + 75.-1 + 71i.o] + \ [71i.i + 7I.i + 7I,-i + 71i,_i] 

(55) 

10 For each n^on-border pixel (represented by a secondary 2-plxei), an 

equation in the form of Equation (55) is obtained. For the border pixels, T = 
l(x) is set. Solving this system, the smoothed image l(x,t) = T Is obtained. 

Optical Flow 

15 ' 

An Indirect application of the heat transfer equation Is the computation 
of optical flow for a 2D Image sequence l(x,t), using the Horn and Schunk [29] 
algorithm. It can be shown that the velocity vector u(x,t) = (u(x,t), v(x,t)) 
satisfies the following constraint arising from variational calculus (for greater 

20 legibility. (x,t) has been dropped) : 

25 

where a Is a weighting factor and Ix, ly and It are the first derivatives of l(x,t) in 
X, y and t, respectively. Let us rewrite Equation (56) in the following vectorial 
form: 
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V/(V/ . u) V^u - Jt V J 

where V^u = (V^u, V^v). Reorganizing the terms of Equation (56), the 
following equation is obtained: 

5 

a^V^n = V/( V/ . u) + VJ. 

(Of) 

Taking C7(x, t) = -VI(VI . u) - It VI as a heat source, Equation (57) 
can be seen , as a steady state heat transfer equation in which the cochain T 
10 corresponds to u(x, t) and Jl=^a^. It can thus be decomposed using the 
method described hereinabove. The cochain T is u(x, t). and the following 
relation is obtained: 

15 

For the same reasons as in the linear diffusion problem, special 
considerations are needed at the borders of the image. Zero velocity is 
assumed at the borders of the image and the system is solved to get the 
velocity field for each point of the image. 

20 

Nonlinear Diffusion 

Linear isotropic diffusion reduces noise but also blurs edges. As the 
scale increases, edges tend to be harder to identify [43]. One possible way of 
25 reducing this effect might be to consider the heat conduction coefficient X as 
a field function dependent on the magnitude of the edges; that is: 



|:/(x,i) = V . (5(|V/(x,t)|^)V/(x,i)) 



(58) 
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Which corresponds to Equation (28) with ;i(x,t) = g(| V l(x,t)|^), T(x,t) = l(x.t), 
p(x) = 1, c = 1 and <T(x,t) = 0 (i.e., unsteady transfer >Anth no source). The 
conduction function g(s) must display the following behavior: in constant 
regions, there should be linear isotropic diffusion (Equation (50)), that is 
5 g(l Vl(x,t)|^) = 1 for I Vl(x,t)|^ = 0, and almost no diffusion when the magnitude 
of the edge is great; that is, g(l Vl(x,t)|^) = 0 for | V l(x,t)l^^oo . Perona and 
Malll< [38] proposied the following functions: 

g(s) = e fc% {k >0) 

The parameter k in these functions is difficult to set because it controls 
15 the threshold of diffusion but also the steepness of the function [35]. An 
advantageous altemative is to use the function: 

20 where k and /control the threshold and the steepness, respectively. Equation 
(58) is then soved for a particular t (the scale) with initial conditions: 

J(x, 0) = /(x) 



10 

and 



25 



where l(x) is the original Image. 
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Let US assume that we have I time steps At = t/l. FIret. the same 

cubical complexes K"' and K' are used as hereinabove and the following 
relations are defined: 



= /C*'x[ifc,tfe+i], tk = kAt, \fke[OJ-i\ 

< 

Secondly, the following assumption is made about the spatial behavior 
of h(x, t); that Is. the approximation function h{x, t) is chosen. For a 3-pixel /e 
as defined In Figure 17. it is assume that H is the mean value over 
[-A/2,A/2] X [-A/2,A/2]. Thus. e = H; that is. using Equation (34). the 
following relation can also be written: 

(59) 

For the sake of simplicity, the same spatial bilinear approximation 
function g (X, t) as hereinabove Is used. The behavior over a time step has to 
be approximated. Some common assumptions about time variation may be 
generalized by proposing [37]: 



A{t) dt = + (1 _ w)A{tk)) At, 0<w<l 



where A(t) is some quantity and w Is a weighting factor [37]. Some values of w 
lead to well-known schemes: 1) w leads to the explicit scheme; that Is. the 
value at tk prevails for the entire time interval except at time tfc+i. 2) w = 1 leads 
to the fully Implicit scheme; that is. the value changes at time t^ from A(tk) to 
Aik*^) and stays there throughout the whole time interval. 3) w = 0.5 leads to 
the semi-implidt or Crank-Nicolson scheme; that Is. there is a linear variation 



10 
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of A(t). It is proposed to use the implicit scheme because for large values of 
. it best emulates long term time behavior for heat [37]. That is for a 3-plxel, 
and for w = 1 : 

g(x) =~[{gI + ^^^^y) ' ^+ (Ol + ^^^x) . j] At, X € 7,. * e [0, At] 

{ 

In order to obtain the local function q (x, t) Equation (25) is applied:' 

q(x,i) = A(x,t)g(x,t) 

wiiere A (x, t) = g(|VI(x.t)|2). As Vl(x,t) is a spatially sampled image where . 
samples are located at the 0-pixels of K" , the local values of X (x, t) are 
approximated. For the sake of simplicity, a bilinear approximation is once 
again used; that is: 

A(x) = a + 6a; + + dxy 

For a 2-pixel of K" , as illustrated in Figure 18. the folllowing relation 
is obtained: 

t 

A(x,t) = Ao,o + ^ ((Ai,o - Ao,o)a; + (Ao,i - Ao,o)y) 4- ^(Ao,o - Ai,o - Ao,i + Ai,i). 

(60) 

Using these assumptions, the same steps as in hereinabove are 
25 followed to find Q as a function of G. For instance, this function for one n-pixel 
as.defined in Figure 16c is: 



15 



20 



Qi = (CL)*G 
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where C, L and G are matrices defined as: 



= 0 



1/24 7/24 1/24 1/24 7/24 1/241 

0 1/24 1/48 0 1/24 1/48 , 
1/48 1/24 0 1/48 1/24 0 J 



Ai,o 



and G 



■I 



'4 



Using Equations (40) and (43), we can express D can be expressed 

as a function of T. For one 3-pixel, of as defined in Figure 19. and the 
following relation is obtained: 



(61) 



where C, L and T are matrices defined as 



1/24 
1/48 
0 

1/48 
-1/12 
0 
0 
0 
0 



1/16 0 

1/2 1/48 

1/16 1/24 
0 0 



1/48 
0 
0 
0 



1/16 1/12 0 0 

0 6/24 0 0 

0 1/12 1/16 0 

1/4 6/24 0 1/48 

-3/8 -7/6 -^/8 --1/12 

5/24 1/4 0 



1/16 1/12 
0 5/24 



1/12 1/16 



0 


0 ■ 


0 


0 


0 


0 


0 


0 


-3/8 


-1/12 


0 


1/48 


1/L6 


0 


1/4 


1/48 


1/16 


1/24. 





























Ao,o 
Ai,o 

- . 


and T ss 



^1 . 



Equation 33 with <5',/^ > = 0 is applied to obtain, for each 3-pixel of 



K' 



7So-(CALA)*TAi = 75 



which defines the system of linear equations. The initial conditions are T° = 
l(x). 
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A Different Hypothesis for Heat Conduction 



In the preceding discussion, A(x) has been approximated with a 
bilinear function, essentially for the sake of simplicity. Nevertheless. It could be 
preferable to use another assumption. Actually, this simple approach does not 
accurately handle abrupt changes in conductivity. For example, two 2-plxels of 
K\ as shown in Figure 20 will be considered. To compute Q, A(x) Is 
apprqximated at tiie borders of the pixels based on the values at their centers. 
Using bilinear approximation, the value of I(x) on the line linking two points is 
declared be the arithmetic mean of the values at these points. For instance, 
given \o and \o -> 1 . the conduction at the border is about 0.5. This 
means that the zero conductivity at one pixel is partly cancelled out by the fact 
that on the pixel beside it. there is a high conductivity coefficient. In non-linear 
gray level diffusion, we are confronted with precisely this kind of abrupt 
change. For example, at step edge pixels, the conduction may need to be very 
low. whereas imnnediately adjacent, it may needs to be almost one. 

A better assumption would thus be to consider I(x) as constant over 
one single 2-pixel of K' [37]. Therefore, on the Lface common to two pixels 
as in Figure 20: 



It can easily be seen that when ;i„.o ->0. then I->0 and when 
\o«Ko' then X^2^o- This means that In both situations, the low 
conductivity would prevail at the boundary common to the two pixels [37]. With 
this assumption, the matrices Q and £^ are modified as follows: 
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ri/4 


1/4 


0 


0 • 


3/2 


-1/4 


-1/4 


0 


1/4 


0 


1/4 


0 


-1/4 


3/2 


0 


-1/4 


-3/2 


-3/2 


-3/2 


-3/2 


-1/4 


0 


3/2 


-1/4 


G 


1/4 


0 


1/4 


0 


-1/4 


-1/4 


3/2 


. 0 


0 


1/4 


1/4 J 



and Ijx = Ao,o 



Ao,^i/(Ao,-»i + Ao,o) 

Ai,o/(Ai,o + Ao,o) 
. Ao,i/(Ao,i +Ao,o) 



Experimental Results 



5 The proposed approach was tested on real and synthetic images in 

the context of linear isotropic diffusion, optical flow and non-linear diffusion. 
The results were compared with another method in each case. 

For linear diffusion, Figure 21a presents our physics-based method 
10 (PB) at three different scales and Figure 21b shows the result by convolution 
for the same scales. In the absence of a quantitative evaluation, it can be said 
that subjectively the results seem to be similar. 

For optical flow, Figures 22a-22c show the first frames of three 
15 sequences: rotating sphere, Hamburg taxi and tree sequences. The results are 
compared with those being obtained using a finite-difference Implementation of 
the Horn and Schunck algorithm (FD) [18, 19]. In these three examples and for 
both the PB and FD methods, the Image derivatives are computed by 
convolution with the appropriate Gaussian derivatives. Both temporal and 
20 spatial scales are set to 1, as is the weighting factor a . Figures 23a and 23b 
shows the flow pattem computed for the sphere sequence. Figures 24a-b and 
25a-b present the flow patterns for the taxi and tree sequences, respectively. 
For the rotating sphere and the taxi sequences, we obtain similar results with 
both methods. For the tree sequence, we also obtain similar results even if the 
25 extreme values seem to be smaller with the PB method than with the FD 
method. This fact is more apparent in Figure 27a and 27b which show 
respectively the results for the PB and FD methods for the tree sequence in 
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which white noise (standard deviation of 10) has been added (see Figure 26). 
Another advantage of the method according to the present invention is that It 
avoids iterations since the algorithm is applied only once. 

5 For nonlinear diffusion, Figures 28a-28c compare the PB method with 

constant hypothesis on X and the FD [38, 17] method for a small window of 
the peppers image with cr = 5. Figure 28b presents the original section with 
white noise added (standard deviation of 10). Figures 28b and 28c show 
respectively the results for the PB and FD methods. One can notice that some 
10 details are better conserved in Figure 28c than in Figure 28b. This fact is 
enlightened in Figures 29a-29d which show a profile of the diagonal line 
starting at the upper right corner and finishing at the lower left corner. 

Figures 30a and 30b present the results for the peppers image with cr 
15 =1.0, 3.0 and 5.0. The results for PB seem a little sharper than the FD 
results. 

Figures 32a and 32b show the results for the Lena image (Figure 31a) 
with an added white noise of standard deviation 10 (Figure 31b) at two 
20 different scales, cr= 4.0 and 8,0. Again, the PB method seems to give 
sharper results at both scales. 

Figures 33a and 33b present details of the Lena results at cr = 8.0. 
Figure 33b seems smoother In constant zones but some details are lost. For 
25 example, compare the eyes, the trim on the hat and the right side of the face. 

Conclusion 

An alternative approach to the PDE-based resolution of the diffusion 
30 problem was described. The proposed approach differs in two significant ways 
from with the classical PDE resolution scheme: 1) the Image Is considered as 
a cubical complex for which algebraic structures such as chains, cochains. 
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boundaries and coboundaries are defined; and 2) the diffusion problem is 
decomposed Into conservative and constitutive basic laws, each of which Is 
represented by cochalns and coboundaries. 

The conservative basic laws are represented without approximation 
while some approximations are required for the constitutive laws. This means 
that unlike traditional PDE resolution, for which many approx.matior,s must be 
made, all approximations are known since they are only needed in the 
representation of the donstltutlve equations. Coboundaries are computed 
usmg fundamental theorems of calculus such as the Green. Stokes and Line 
Integral theorems. Unlike Iterative numerical analysis algorithms that do not 
allow the explanation of intemiedlate results, the use of basic laws allows the 
Physical explanation Of all steps and intem^edlate results of the algorithm 
Moreover, since there is no Iteration In the resolution process, there Is no 
problem about the convergence of the numerical analysis scheme 
Furthemiore. the use of cubical complexes provides algorithms that can 
operate In any dimension. It has the significant advantage of avoiding the 
potentially difficult task of extending the algorithm to higher dlr^ensions 
Cochalns and coboundaries allow the use of both global and local quantities" 
Integrals or discrete summations over fields are used to compute global 
quantities. This allows the reduc^on of noise by performing a smoothing 
operation, as opposed to differentiation, which enhances high frequencies. 

In computer vision and Image processing, several problems can be 
modeled as diffusion problems. The proposed approach has been validated on 
smoothing by linear and nonlinear diffusion and on the computation of optical 
flow. The results obtained conflmi the effectiveness of this approach 
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PRACTICAL EXAMPLE #2: 



A PHYSICS-BASED MODEL FOR ACTIVE 
CONTOURS 
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A new active contours model is presented. It is based upon a 
deconriposition or the linear elasticity problem Into basic physical laws. As 
opposite with other physics-based active contours model which solve the 
partial differential equation arising from the physical laws by some purely 
5 numerical techniques, exact global values are used and approximations made 
only when they are needed. Moreover, these approximations can be made 
wisely assuming some l<nowledge about the problem and the domain. The 
deformations computed with the present approach have a physical 
interpretation. In addition, the deformed curves have some interesting physical 
10 properties such as the ability to recover their original shape when the extemal 
forces are removed. The physical laws are encoded using the computational 
algebraic topology based image model described herein. The resultant 
numerical scheme is then straightforward. The image model allows our 
algorithm to perfomn with either 2D or 3D problems. 

15 

Introduction 

These last years, active contours and active surfaces have been 
widely studied since the introduction of active contours by Kass et a! [59]. They 
20 have been used in image segmentation [62], tracl<ing [68], automatic 
correction and updating of road databases [46], etc. 

To solve these problems, many dilTerent approaches have been 
proposed (see [57, 63]) in particular physical models derived. from equations of 

25 continuum mechanics. Mass-springs models are physical models which use a 
discrete representation of the objects. Objects are modeled as a lattice with 
point masses linked together by springs [57]. Infomnation is thus only available 
at a finite number of points [63]. These methods offer only a rough 
approximation of the phenomena happening in a body [57]. Moreover, the 

30 determination of spring constants reflecting the material properties may be a 
very fastidious work. However, they offer real-time performances and allows 
for parallel computations. 
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Other physical models based upon the minimization of an energy 
functional which takes into account an internal regularizing force and an 
external force applied on the data are also often used. Some of them consider 
the deformable bodies as continuous objects by approximating their 
continuous behavior with methods such as the finite element method (FEM). 
FEM are closer to the physics than mass-springs models but their 
computational requirements make them difficult to be applied in real-time 
systems without preprocessing steps [57]. Finite difference methods (FDM) 
are also used to discretize the objects. They usually offer better performance 
than FEM but they require the computation of fourth order derivatives which 
make them sensitive to noise [63], 

For a given curve S, the application of the FEM and FDM methods 
leads to a discrete stationary system of equations of the form: 



where K is a matrix which encodes the regularizing constraints on S and f(S) 
represents the data potential. However, some problems such as animation in 
graphics applications require to take into account a dynamic evolution of the 
curve [57]. In these case, inertial body forces and damping forces may also by 
considered by controlling the deformations through a Newtonian law of motion: 



Q2C QC 



or by a Lagrangian evolution 



(63) 
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Where M and D are respectively matrices which represent the mass model and 
the background damping. Equations (62) and (63) are solved using various 
numerical schemes [70. 64] assuming an Initial curve So close to the solution 
which evolves until the inertial tenns go to zero. 

Over the last years, a lot of different methods have been introduced to 
compute the matrices IVI. D and K but as pointed out by Montagnat et al [63], 
these methods have a major drawback: the con-esponding system 
deformations do npt have any physical Interpretations. 

A new model which includes a physical interpretation of the 
deformations is described hereinbelow. The model is similar to a mass-springs 
model but It includes both the efficiency of the mass-springs models and the 
accuracy of the physical modeling of the FEM by providing a systematic 
method for specifying springs constants reflecting the properties of the 
materials. 

To achieve it, we propose to use directly the basic laws of physics 
which lead to the partial differential equations (62) and (63). These equations 
are Indeed obtained by considering and mixing together some basic laws of 
physics into a global conservation law and considering its local counterpart 
[72]. This approach is not always well suited for problems such as computer 
vision in which the continuous domain must be subdivided into many sub- 
domains for which there are often only one information available. The use of 
this Infomriatlon as a global value over each sub-domain allow to directly use 
the global conservation law which can lead to an algorithm less sensitive to 



noise. 



To encode these global values over points, surfaces, volumes, etc 
arising from some physical laws, we use the computational algebraic topology 
based image model described hereinabove. 
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The approach according to this illustrative embodiment of the present 
invention has several advantages. 1) Since the linear elasticity problem is well- 
known in continuum mechanics, the modeling according to the illustrative 
5 embodiment of the present invention can be made wisely in order to provide 
some good physical interpretation of the whole deformation process and of its 
intermediate steps. This allows an easier determination of the parameters, 
used in the process since they have a physical meaning; 2) The determination 
of the springs constants in order to reflect the material properties is 

10 straightforward; 3) The objects in the image (e.g. curves, surfaces) are 
modeled as entities having their own physical properties such as elasticity and 
rigidity. They have the property of recovering their original state when the 
forces applied on them are removed; 4) Both smooth results and results 
having high curvature points can be obtained; 5) The complexity of the 

15 algorithm Is minimal and allows for real-time simulation without any extra 
preprocessing steps [51]; 6) The image model allows our algorithm to perform 
with either 2D or 3D problems. 

Physical Modeling 

20 

One of the objectives of this illustrative embodiment of the present 
invention Is to model the objects in an image (e.g. curves, surfaces, etc) as 
entities having their own physical properties such as elasticity and rigidity. As a 
consequence, these objects need to satisfy the laws and principles of the 
25 continuum mechanics. For Instance, a body subjected to forces must move or 
deform according to the universal laws of physics. 

These principles and laws to which all bodies must obey will now be 
presented. We first Introduce the concepts of stress and strain which are 
30 required in the statement of the governing equations for deformable bodies- 
Then, the physical laws related to the linear elasticity problenn will be 
presented. 
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The elasticity tlieory has been widely studied by engineers and 
scientists and is the main subject of many books. The present specification 
only presents the concepts of this theory which are relevant to our application. 
5 The concepts presented here are well known and may be found in many 
continuum mechanics books such as [65, 50]. 

Forces. Stresses and Strains 

A material body in a 3-D space is always subjected to forces. These 
forces may come from an external agent (external forces) or issue from the 
object itself (internal forces). When the externa! forces are greater than the 
internal forces, then the body can undergo deformations (strains) or be 
accelerated. This deformation can induce internal forces (stresses) if the 
material is elastic. The concepts of force, stress and strain and the relation 
between strain and stress is exposed hereinbelow. 

Forces acting on a body 

Two basic types of forces act on a body. First, there are the interatomic 
forces which hold the body's particles together at some configuration. These 
forces, called internal forces, can either attempt to separate or bring the 
particles closer according to the fact that the body undergoes a contraction or 
an extension [65]. They act in response to a force applied by some other 
25 agent. Assuming the equilibrium of the body and using Newton's law of 
reactfon, they must be equal In magnitude to the forces applied to the body but 
in opposite direction [66]. 

On the other hand, there are the forces applied by an external agent. 
30 called external forces. Two types of these forces are generally applied on a 
body. First, there are forces such as gravity and inertia called the body forces, 
which act on all volume elements. These forces, noted bi (forces per unit of 



10 



15 
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mass in a direction Xi), are distributed In every part of the body. Secondly, 
there are forces which act on and are distributed over a surface element such 
as the contact forces between solid elements [49]. These forces, noted fi 
(forces per unit of area in a direction X|, are called the surface forces. The 
5 surface element may be inside the body or a part of a bounding surface [60]. A 
body of arbitrary size, shape and material subjected to surface and body 
forces Is shown in Figure 34. 

External forces applied on a body must be transmitted to it, A rigid body 
10 can then undergo either a spatial shift, a rotation of both of them. In the case 
of a non-rigid body, it can also go through a deformation or a distortion in 
which case internal forces will be developed to counterbalance the external 
forces. If the internal and external forces are balanced, we say that the body is 
In static equilibrium. Otherwise the body can be accelerated which would give 
15 rise to inertia forces [58]. Using d'Alembert's principle, these forces may be 
included as part of the body forces [48] such that the equilibrium equations can 
be satisfied. If the body gets deformed then the deformation can either be 
elastic or not and is subjected to the material properties of the body such as its 
elasticity and its rigidity. If the internal forces induced by the material 
20 properties of a body are too weak to counterbalance the extemal forces, then 
the body can be pemianently defonned [52]. 

We assume herein the material to be isotropic with respect to some 
mechanical properties. We then suppose that the material properties are the 
25 same in all directions for a given point [60]. We also consider an 
homogeneous material which means that Its properties are identical at all 
locations. 

The Concept of Stress at a Point 

30 

Let us consider an isotropic and homogeneous body B. Let us assume 
that B is subjected to arbitrary surface and body forces such that B Is in static 
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equilibrium. Let P be a Interior point of B and S be a plane surface passing 
through P. S will be referred to as the cutting plane and is defined by the unit 
nomrial vector n=(ni, na, na)"^. Then S partitions B into two sections I and II as 
shown in Figure 35. 

5 

Let us assume that AS \s a small element of area of the cutting plane 
sun-ounding P (see Figure 36). 

Since the body is in static equilibrium then the force system acting on 
10 each part I and 11 taken alone must also be in equilibrium. This generally 
requires that some internal forces are transmitted by part I to part II. These 
forces are not necessarily distributed uniformly on every part of the cutting 
plane. Thus they may vary in magnitude and direction over It. We generally 
want to determine precisely that force distribution at every point of AS . The 
15 temn stress is used to define the Intensity and the direction of the internal force 
A/" acting at point P. Using the Cauchy stress principle [60] we define the 
stress vector (or traction vector or traction forces) t" at P as: 



t"= lixn 



20 



assuming that P remains an Interior point of AS as Its area reduces to zero. 



Let us mention that t" Is not necessarily In the direction of the normal 
vector n at P. However. It may be decomposed into a component 

25 perpendicular to the cutting plane, called the normal stress, and a component 
parallel to It. called the shear stress. The nomnal stress attempts to separate 
(bring closer) the material particles after a compression (an extension) of an 
elastic body when it tries to recover its original state. On the other hand, the 
shear stress acts parallel to the cutting plane and tends to slide adjacent 

30 planes with respect to eadi other (see Figure 37a-37d). 
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It Should be first noticed that the stress vector e is defined with respect 
to the cutting plane's unit normal vector n. Since there are Infinitely many 
cutting planes going through P there are also as many stress vectors defined 
at P. Juvinall [58] defines the state of stress at a point P as a complete 
description of the stress magnitude and direction for ail possible cutting plane 
passing througfi P. Fortunately, this description can be fully obtained by 
considering any tlnree mutually perpendicular planes passing at P [58]. For the 
sal<e of simplicity, we usually use the three axes defined by the three 
canonical vectors Xi, Xa and X3. 

Let us define <r. as the stress component in the direction of when the 
nonmai vector is parallel to the axe defined by xi. If I = j then represents a 
nomiai stress. Othenwise, cr^ is a shear stress. With these conventions, the 
component t," in the direction of x, of the traction force t" depends on the 
nomial stress o-„ , the shear stresses <r^,and o-„ and the nonnai vector n such 
that (see Figure 38): 



if = cfuni + (72in2 + (J^ns 

3 

^'=^ (64) 

20 

Equation (64) is known as the Cauchy stress fonnula. 

Since each of the three coordinates axes involves six stress 
components, there is a total of nine stress components. However the 
25 equilibrium of monnents at P [49] gives that only six of these are independent 
that is ffg =<Tj,for all i. J = 1, 2. 3. This means that the state of stress at a point 
is fully detenmined by o-„, <7^, a^,, a,^, cr^. 
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77ie Concept of Strain at a Point 

Any non-rigid body goes through defonnations and distortions when 
subjected to forces. The body can either extend or contract (deformation) or 
5 have a geometric modification of its shape (distortion). Figures 39a and 39b 
present these concepts. 

The term strain refers to the direction and intensity of the deformation 
at any given point with respect to a specific plane passing through that point 
10 [58]. As for stress, the strain is defined according to a specific cutting plane. 
The state of strain Is defined by Juvinall [58] as the complete definition of the 
magnitude and direction of the deformation at a given point with respect to a 
all cutting planes passing through that point. 

15 As for the state of stress, the description can be obtained by 

considering any three mutually perpendicular planes passing at P. One can 
therefore see a great similarity between stress and strain. However, there is a 
major difference between them: strains are generally some directly 
measurable quantities while stresses are not. Fortunately, stresses can be 

20 computed from strains (and vice versa) using a constitutive equation. 

As for stress, two types of strains can be defined. First, there are the 
strains which result of a change In the dimensions of the body (deformation). 
Let B be the same body defined hereinabove, AB be a small element of B of 

25 length Ax, in the direction of xi and AB' be the deformation of AB such that 
AUf is the change in lenglit of AB after the application of a force in the 
direction of Xi (see Figure 40). The normal strain at P In the Xi direction with 
respect to a cutting plane having xi as nomrial vector is the unit defonnation of 
a line element [65] in the direction of Xj. It Is formally defined as: 

30 ■ 
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6u = lim 



Aasi-^O Axi dXi 



The normal strain is clearly the unit change in length per unit original 
length for the element in the direction of Xj, Since it is the ratio of two units of 
5 length, it is dimensionless even if it is sometimes expressed as units of length 
per unit of length such as inches per inch. 

Let us now suppose that there are two perpendicular lines PB and PA 
of length Axj and Ax^^ respectively In the direction of Xj and Xk (see Figure 
10 41). 

Let us assume that after a distortion points A and B move respectively 
to A' and B' while P remains fixed. The lines PA and PB have been rotated of 
angles djj^ and such that: 



15 



20 



25 



^ = tan{ejk) and ^ = tBJi{0kj) 

where Au^. and Au^ are respectively the displacements of B and A in the xj and 
Xk directions. 



If it is assumed that only small distortions occur, then we can 
approximate both tangents by their angles. Thus: 



or, taking their infinitesimal analogous: 



10 



15 
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^ ^ (65) 

The shear strain at P with respect to the cutting plane having x, as 
normal vector is the angle in radians through which two orthogonal lines in the 
undistorted body are rotated by a distortion [56]. That is = Oj^ + . The two 
subscripts in have a similar meaning as for stress. For instance. r« is the 
strain acting on two adjacent planes perpendicular to the x, axis and sliding 
them relatlve to each other in the Xk direction. 



Let us recall that these definitions have been made under the 
assumption that only very small displacements occur in the body. The normal 
and shear strains are supposed small compared to unity {50]. If this constraint 
is relaxed in order to Include large defomiations then the system to solve for 
the computation of the forces, the stresses, the strains or the displacements 
becomes non-linear and then harder to solve. This is sometimes necessary in 
some problems where large deformations can occur such as for thin flexible 
bodies [50] of for the modelization of human tissue [53]. However, if we restrict 
ourselves to small defomiations. then the approximations made to define the 
shear strains do not induce too many errors and are widely accepted in the 
20 dasslcal theory of elasticity [69, 65. 49. 50]. 

Finally. Equation (65) dearly shows that the shear strain is also a 
dimensionless quantity since it is the ratio of two units of length. Defining for i ^ 
j: 

25 



SiJ = piJ (*,J = 1,2,3) 
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25 



the followin strain-displacement relation (or tlie Idnematical relatlonsiiip [69]) is 



As for stresses, there is a total of nine strain components at each point 
of the body (three per mutually perpendicular cutting planes) but by symmetry 
they can be reduced to six. that is /jt = r^- for all j. k = 1 , 2, 3 with j * k. The 
state of strain at any point can then be described by , , s^^ , s,^ , e,^ and 



Relations Between Forces, Stresses and Strains 

As mentioned hereinabove, strains are measurable quantities while 
stress are not even if both can be computed from the other. This is due to the 
fact that some knowledge about the material of a body Is necessary to 
measure the stresses from strains and vice versa. For instance, a steel beam 
and a rubber beam Induce different Intemal forces when bent equally. 

This gap between strains and stresses will now be filled by stating the 
physical laws relating them to each other. This hole between strains and 
stresses needs to be filled by a constitutive equation (or material law), which 
refiect the intemal constitution of the materials. The material law for the linear 
elasticity problem Is known as the Hooke's law. 

Before stating the law. It should be first reminded that a material has an 
elastic behavior when it satisfies the two following conditions: 

1 . The stresses depend only on the strains. 

2. Its properties allow a body to recover its original shape when the external 
forces applied on the body are removed [60]. 



obtained: 




(66) 



5 



s. 
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If theses conditions are not satisfied, a body is said to have an Inelastic 
behavior. Any body may be seen as having an elastic behavior as long as It Is 
not defomned beyond a limit value. This value Is called the elastic limit [52, 49] 
5 and Is usually defined as the maximum value of stress that a body can 
undergo without undergoing a peananent deformation. 

In addition. If the stress is a linear function of the strain, an elastic 
material has a linear elastic behavior. In what follows, it is assumed that the 
10 material has a linear elastic behavior. 

As mentioned in [60] and [49], in many situations the problem of 
elasticity can be considered as a 2D problem. This particular case Is known as 
plane elasticity. Two basic types of problems compose the plane elasticity. 
The problems in which the stress components in one direction for a body are 
all zero are referred to as plane stress problems. On the other hand. If all the 
strain components In one direction for a body are zero then the state of strain 
for that body is refen-ed to as plane strain (see [65, 60, 58]). 



15 



20 



25 



The present problem is considered as a plane strain problem. This 
distinction is important since the constitutive equation slightly differs according 
to the fact that a plane stress or a plane strain problem Is considered. 



The Hooke's Law 



When a rubber ball Is compressed its diameter In the directions 
perpendicular to the applied force gets larger. A similar phenomenon occurs 
when a rubber band Is extended and Its cross section gets smaller\ In fact, 
these changes In dimensions happen In all materials even If they can't always 
30 be noticed by a naked eye [52], 



* Example taken in [52] 
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When a stress is acting on an isotropic and homogeneous body in only 
one direction (uniaxial stress), one can show that the transverse strain (the 
strain in a perpendicular direction) is directly proportional to the strain e 
5 induced by the stress: 

= —ve 

The ratio: 

10 




is called the Poisson ratio, it is supposed constant when the stress is below 
the elastic limit [66]. 

The linear relationship between a uniaxial stress cr^ in a direction xi 
and the con-esponding strain is l^nown as the Hooke's law [52. 58] and is 
written as: 



where E is the Young's modulus of elasticity. Of course, the Hool<e's law is 
valid if the stress is not beyond the elastic limit of the material. 

As pointed out by Boresi [49], the stresses at any point depend on all 
the strains in the neighborhood of that point. Thus the total deformation in the 
X| direction depends not only on the stress in that direction but also of the 
defonnations in the other two perpendicular directions. For instance the nornial 



wo 2004/015631 



86 



PCT/CA2003/001214 



15 



20 



strain does not only depend on ^(Hooka's law) but also of the transverse 

E 

strains s-^- and e^^ such that the total deformation in the direction of X| is: 



E E E 
E 



k« - {cTjj + afcfe)] {% T^j^k) 

(67) 



Equation (67) is the normal strain-stress relation. For isotropic 
10 materials, it can be shown that the normal strains are not Influenced by the 
shear stresses [52]. Consequently, the shear stresses only induce shear 
strains and they are related by the relation: 

^ (68) 



where G = ^^^^^^ is called the modulus of rigidity. Equations (67) and (68) 

are known as the Generalized Hooke's law for linear elastic isotropic materials 
[52]. These equations can be inverted in order to express the stresses as 
functions of the strains: 



E 

=^ (1 + u){l - 2v) - ''^^^ + ""^^^^ 

(69) 



E 

^ (70) 
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Relation Between Forces and Stresses 

It is well known that the conservation (balance, equilibrium) laws 
5 constitute an Important class of equations in continuum mechanics. They 
relate the change in total amount of a physical quantity inside a body with the 
amount of this quantity, which flows through its boundary. These laws must be 
satisfied for every continuous materials. Local differential equations are usually 
used to express these laws. In what follows, the linear momentum, which is 
10 relevant for the linear elasticity problem, is presented. 

To every material body B Is associated a measure of its inertia called 
the mass. This. measure may vary in space and time inside a body. Let V be 
the volume of B, S its bounding surface and Am be the mass of a small 
15 amount of volume AV . The mass density is given by: 

p = p(x,<)= lim 

Let us assume that distributed body forces /?t&. and tractions forces t" 
20 are applied to S (see Figure 42). Let also assume that B is moving under the 
velocity field vi = V|(x,t).. The quantity: 

Pi(t)= jjj puidV ' ' 

V 

25 is called the linear momentum of B. The principle of linear momentum [49]. 
states that the resultant force acting on a body is equal to the time rate of 
change of the linear momentum. Thus: 
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V 8 V 

> s.,- ' 

Forces acting on the body j 

Recalling Equation (64): 

3 

where n = (ni, n2, x\zf Is the unit nomnal vector to the surface and 
<T = (cri.,cr2,,cr3.)^and using Gauss's divergence theorem, the following 
relation is obtained: 



/// rfv = /// <iK = /// ( V . + A) rfv 



^ (72) 



Since V is arbitrary tlien tlie integral sign can be retrieved leading to \he 
local equations of motion: 



(73) 



The global equilibrium equations can be obtained assuming a zero 
velocity field in Equation (72): 

20 
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fff'^'O^ rf^ + fff P^i 0 (i = l,2,3) 

V V 

V ' ^ V ' 

iQCeraal foices External forces 



and their local counterparts are: 

V-a + pbi=0 (2 = 1,2,3) 



(74) 



(75) 



Let us now summarize the decomposition of the problem. The relation 
between the global displacement U of a body and the corresponding strains 

10 (see Figure 43a) has been introduced hereinabove. The constitutive equation 
relating the strains and the stresses (see Figure 43b) has also been 
presented. Finally, how the stresses are related to forces using the linear 
momentum principle (see Figure 43c) has been described. A general scheme 
similar to the one presented by Tonti [72] may then be introduced to 

15 summarize how the internal reaction forces of a body are related to the global 
displacements of that body (Figure 44). 

Discrete Representation of Images 

20 Some algebraic tools used to model the image will now be recalled 

from the above description. An image is composed of two distinctive parts: the 
Image support (pixels) and some field quantity associated to each pixel. This 
quantity can be scalar (e.g. gray level), vectorial (e.g. color, multispectral, 
optical flow) or tensorial (e.g. Hessian). The Image support Is modelled in 

25 terms of cubical complexes, chains and boundaries. With these concepts, it is 
possible to give a formal description of an image support of any dimension. 
For quantities, the concept of cochains which are representations of fields over 
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a cubical complex is Introduced. For the use of these concepts In image 
processing, see [45]. 

An image Is a complex of unit cubes usually called pixels. A pixel 
5 y c: 91" is a product : 

7 = /l X /2 X . . . X 

where Ij is either a singleton or an inten^al of unit length with integer end 
10 points. Then Ij is either the singleton {k} and is said to be a degenerate 
interval, or the closed interval [k, k+1] for some k^Z. The number 
^e{0,l,...,n}of non-degenerate intervals is by definition the dimension of 
which is called a q-plxel. Figures 45a-45b illustrate three elementary pixels in 

15 

For q>l, let J = {^o,^„...,^^,}be the ordered subset {1, 2 n} of 

indices for which I^^ =[aj,bj\ is non-degenerate. Define: 

AkjCT = X . . . Ikj^x X {o^} X /jy+i X . . . X Jn 

20 

and 

25 The and tine B^^ are called the (q-1)-faces of <r . One can define the 

(q-2)-faces. .... down to the 0-faces of <t the same way. The faces of y 
different from / itself are called its proper faces. 
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By definition, a natural orientation of the cube is assumed for eacti 
pixel. Suppose that y denotes a particular positive oriented q-pixel. It is 
natural to denote the same pixel with opposite orientation by -y . Examples of 
orientations are given in Figures 45a-45c. A cubical complex In SH" is a finite 
collection K of q-pixels such that every face of any pixel of the image support 
called K is also a pixel in K and the intersection of any two pixels of K is either 
empty or a face of each of them. For example, traditional 2D image models 
was only considering pixel as a 2D squared element. Definitions presented 
before allows us to consider 2-pixels (squared elements). Lplxels (line 
elements) and 0-pixels (punctual elements) simultaneously. 

In order to write the image support in algebraic form, the concept of 
chains is introduced. Any set of oriented q-pixels of a cubical complex can be 
written in algebraic form by attributing them the coefficient 0.1 or -1 . if they are 
not in the set or if they should or not be taken with positive orientation, 
respectively. In order to represent weighted domains, arbitrary Integer 
multiplicity for each q-pixel is allowed. 

Given a topological space X c 91" in terms of a cubical complex, a free 
20 abeiian group Cc(X) generated by all the q-pixels is obtained. The elements of 
this group are called q-chains and they are formal linear combinations of q- 
pixels [45]. A fomial expression for a q-chain c, is c, = Y,^^A,r, where ^,^z, 



15 



25 



The last step needed for the description of the image plan is the 
introduction of the concept of boundary of a chain. Given a q-pixel y, its 
boundary dy is defined as the (q-l)-chaln con-esponding to the alternating sum 
of its (q-l)-faces. The sum is taken according to the orientation of the (q-1)- 
faces with respect to the orientation of the q-pixel. It is said that a (q-1 )-face of 
y is relatively positively oriented with respect to the orientation of y if its 
30 orientation is compatible with the orientation of . By linearity, the extension of 
the definition of boundary to artjitrary q-chains is expedient. For instance, in 
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Figures 45b and 45c, the boundary of the 1 -pixel a Is X2 - Xi and the boundary 
of the 2-plxel A Is a + fe-c- Jwhereby a and b are positively oriented with 
respect to orientation of A but c and d are negatively oriented with respect to 
orientation of A. Let us notice that the boundary of a 1 -pixel Is always the 
5 difference of Its boundary points. The boundary can be defined recursively. 
Suppose a (q-l)-chain and a q-chain defined as =/^^, x[a,6], the 
boundary of can be recursively written as: 

dj^ = ^7,-1 X [a, 6] + (-l)«'-^)(7,_i X {6} - 7^-1 x {a}) 

(76) 

10 

In order to model the pixels quantity over the innage plane, an 
application F has to be found to associate a global quantity with all q-pixels y 
of a cubical complex. This application is denoted <F,y>. This quantity may 
be any mathematical entities such as scalers, vectors, etc. For two adjacent q- 
15 pixels y, and y^ , F must satisfy < F,\y, + X^y^ >^X,< F,y^>+2^ < F,y^ > , 
which means that the sum of the quantity over each pixel is equal to the 
quantity over the two pixels. The resulting transfonmation F:C^(X) ~>5R is 

called a q-cochain and is used as a representation of a quantity over the 
cubical complex. 

20 

An operator is finally needed to associate a global quantity to tiie (q+1)- 
pixels according to the global quantities given on their q-faces. Given a q- 
cochain F, an operator d , called the coboundary operator, is used to transfonn 
F into a (q+1)-cochain dF such that: 

25 

<6J^,J>=<J^,dj> 

for all (q+1)-chains /. The coboundary is defined as the signed sunn of the 
physical quantities associated with the q-faces of /. The sum is tal<en 
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according to the relative orientation of the q-faces of the (q+1)-plxels of y with 
rasped to their orientation. Figure 46 presents an example of the coboundary 
operation for a 2-pixel. 

5 Representation of the Equilibrium Equation 

The basic laws of Figure 44 with concepts of algebraic topology in 
order to get a generic algorithm for solving the equilibrium Equation (74) will 
now be modeled. 

10 

The algorithm is resumed as follows: 1) The image support is firstly 
subdivided into cubical complexes; 2) Global quantities are computed over 
pixels of various dimensions via cochalns according to basic laws; 3) The 
constitutive equations 69 and 70 are expressed as a linear transformation 
15 between two cochalns. 

Ttie Relative Displacement 

Let B be a body In a 3D space and K' be a 3-complex representing the 
20 subdivided spatial support of B. Let us consider a 0-cochaln u and a 1 -cochain 
D such that D is the cot>oundary of u : 

X>:Ci(X?) -5- M3 
25 Figures 47a and 47b present some examples of u and D for a 3-pixel of 



30 



The computational rules for both cochalns u and D will now be 
spedfied. Recalling the strain-displacement relation (Equation (66)): 
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1 duj dui 



(79) 



we have an application : 

U ^'(U) = (eii 5 £22j fi33) ^12} Sl3> ^23)^ 

Omitting the shear strain cjomponents as In [67], the following relation 
may be defined: 



Using the global form of Equation (80) over a 1 -pixel such that 
dy^^ =x^-x^, the following relation Is obtained: 



where t/;K/> is an infinitesimal part of the domain Since Vu is a 
conservative field, applied is the line integral theorem [55, 71] which states that 
for a conservative field F(x) = Vf(x) and for two points A and B in an opened 
connected region containing F(x) the integral of the tangential part of F(x) 
along a curve R joining A and B is independent of the path: 



(80) 




^ F(x)di2 = /(5)-/(A) 
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From Equation (81). the following relation is tinen obtained:' 



(82) 

5 On tlie other hand, applying the cochain D to the 1 -pixel 'eacls to : 

< ^, 7d >=< djD >= U(X^ - X#) = U{x^) - U(x^) 

(83) 

10 which is the same fomi as Equation (82). U(x) = U(x) is then defined. 
Consequently, the location of the displacement vector U must correspond to 
the 0-pixels of K' . The previous definitions are extended by linearity to the 1- 
chains of jfiT^ . , ' 

1 5 The Force-Stress Relation 

Let us consider another 3-complex K' also representing the subdivided 
spatial support of the body B. Let us also consider a 3-cochain f and a 2- 
cochain S such that F is the coboundary of S: 

20 

7 ^ <:F,j>=<dS,j>=<S,dj> 
Figures 48a and 48b present some examples of f and S for a 3-pixel of 

25 



Let be a 3-pixel of ^*and /s be a 2-chaln over K' such that 
= ^/f • Let us assume that the 2-faces y^j of rp ^r® relatively positively 
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oriented with respect to the orientation of Yf • The definition of the coboundary 
leads to: 

'^•^ (85) 

5 

Again, the computational rules associated with F and S are determined. 
Setting a zero velocity field In Equation (71) leads to: 



V s 
10 

where a-, =(o-,„a-2„cT3,), To fulfill Equation (85), the following relation is 
defined: 

<^<,Ti^>= fJJ pbidV 

y (86) 

15 

and 

<Si,-fs^ >= jj-^i • (i= 1,2,3) 

S (87) 

20 where 



and 
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5 TTie Stress-Strain Relation 

Exact global versions of Equations (66) on K" and (74) on haying 
been presented by the means of Equations (83). (86) and (87). In order to 
complete the scheme of Figure 44, representation of the Hooke's law 
10 (Equations 69 and 70) which links the local values of ^ (x) and o-(x) Is 
needed. Since Equations (69) and (70) are constitutive equations, a 
topological expression thereof cannot be provided. Instead, they are 
expressed as linear transformations between the cochains D and S; 

15 

To find this transfomiation, Equation (87) Is recalled: 

< Si,ysi >= JJ^ "(Ti ' ndS (i = 1,2,3) 

20 

which links the cochain S with the strains e using the generalized Hooke's 
law. Unfortunately, the strains are only known at a finite number of points and- 
are then approxinnated over the whole domain S with an approxinnation 
function ^(x). i'(x) Is chosen such that for each 1-face y^, of a 2-pixei / of 
25 K'' , the following relation Is obtained: 



I e{x) • dR=< V^'Yd > 

ID 



(88) 
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5 



where dR is an infinitesimal part of the domain represented by y^^ . It should be 
pointed out that that only approximation of the normal components of s is 
needed. In fact, given: 



where tr(x) is the approximated displacement vector over y and if f^(x) is 

chosen to satisfy Equation (88), then the vector U (x) is fully determined. Then 
the shear components of ?(x) can be computed by appropriately 
10 differentiating the components of £7(x). Using this remark and applying the 
generalized Hooke's law to s (x) satisfying Equation (88), the following relation 
is obtained: ' 

E 

^ E ^ 

at all point of y . Equation (87) Is then replaced by: 

<Si,ysj>=^ff -3=<(x) . iid5 = Ai(e) (i=l,2,3) 
20 '''' S (89) 

which depends on the choice of the approximation function six) and of the 
relative position of K' with respect to iiT' . 



25 Summary of the Algorithm 
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The algorithm used to find an expression of the internal forces 
according to the displacements of a body is now summarized. The input data 
for this algorithm are the cochain U and the material properties of the body 
(the values of E and v ). 

1 . Choice of the location of with respect to K' . 

2. Computation of the cochain D. 

3. Computation of the cochain S 

(a) Choice of the approximation function e (x). 

(b) Application of Equation (89) to express S as a function of the 
displacement components. 

4. Computation of the force by applying Equation (85). 

Applications 

I 

I 

Active Contours 

The above described approach is applied to a 2D active contour model 
based upon a Lagrangian evolution of the curve S: 

dS 

(90) 

where K is the matrix which contains the regularization forces of the curve. 

This dynamic system is discretized in time using a finite difference 
scheme. For a given time step At. the time derivative can be approximated 
by: 

dt At 
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The curve deformation is governed by each vertex displacements 
compared to their neighbors until the equilibrium between the inertia forces, 
the image forces and the Internal forces. Equation 90 is solved using an 
explicit scheme: 

5 

Assuming that the initial curve So was in an equilibrium state jand that 
the initial body forces Fo = KSo are constant during the deformation process, 
10 these forces can be added to the extemal forces Fext leading to a modified 
version of Equation (91): 

= 5i + Ai(F«,-(Ki9t"K5o)) 

= 5t + At(F«, ^ KU) ' (92) 

15 where U is the displacement vector of the curve S. 

The image subdivision process is similar to the one presented in [57]. 
Here, it is desired to solve Equation (92) for local U(x) located at the center of 
each pixel and known initial curve So closed to the solution. Following the 

20 steps presented hereinabove, the two dimensional cubical complexes and 
K' are first positioned. As mentioned, is placed in order to have its 0- 
pixels corresponding to the center of the image pixels. K' is positioned in 
such a way that its 2-pixels coincide with the image pixels. Thus, the 2-pixels 
of are rectangular and symmetrically staggered with the 1 -pixels of and 

25 each 1-pixel of intersects orthogonally in the middle of a 1 -pixel of JST". 
Mattiussi [61] showed that this way of positioning AT^ allows the use of lower 
order approximation polynomials without losing accuracy. Figure 49 shows the 
two complexes positions for a 5x5 - image. 
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In order to solve Equation (92), global values F are needed over each 
pixel of K' . Since these values are generally known, there Is no need to try to 
express them in a topological way. In the exannples, the gradient field of the 
5 bright line plausibility Image obtained using a line detector proposed In [54] 
has been used. It was assumed that the gradient provides global values valid 
over the whole pixel. Thus F = VL Is set wliere L Is the line plausibility Image, 
VL = g'^*\. and g'^ is the Gaussian derivative at scale <t. An approximation 
function e (x)=(f„ (x), s^^ (x)/is also chosen. For simplicity, It is assumed that 
10 ^ (x) = VC7 (x) arises from a bilinear approximation: 



U(x) = (t4(x), i72(x))^ = a + ba;i + cxi + dxi^a 



Thus: 



15 



e22(x) = c + daji 



20 



Since (x) and e^^{yC) has to satisfy Equation (88) for all 1-faces y^, 
of any 2-plxel y of K" as in Figure 50, the following relations hold: 




25 
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from which the following reiation is obtained: 

= i [(p. + ^x.) . (d. h- ^..)] = vuw ^^^^ 

From Equation 93 and tlie definition of tlie nonnal strains, it is 
straightforward that: 

U(X) A + 1 (2?iXi + V^X2) + ^ (^3 - A + 2>2 ~ V4) Xj^X2 

(94) 

where = c7(0) . Equations (83) and (94) lead to: 

ij(x) = fj(a:i,X2) = U(0,0) + |-(U(A,0)-U(0,0))a;i + ^(U(0,A)-U(0,0)) 
+X2 (U(0, 0) + U(A, A) - U(0, A) - U(A, 0)) 3:10^2 

(95) 

from which the values of a, (x) = ( 5^,, (x), o^^. can be deduced. 

The last step is the computation of the internal forces F for each 2-pixel 
of K'. With K" ahdK' positioned as mentioned before, each 2-pixel /j, of 
K' intersects four 2-pixels y^. y^, ;.^and of K" . That is. four 
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approximation functions a,!^ , , af and Srf corresponding to the four 
Intersecting 2-pixels of (see Figure 51) have to be considered. 

The value of the cochain S are found over the four 1-face of yby the 
appropriate integration: 







• ~^ dx2 






' — J dxi 






' ~f dxx 


^ = i'-^(-f-)• 


-l*dx, + /%f(-|,.,) 


• — 1 aX2 



Equation (85) leads to: 



By substituting Equations (96) and (95) in Equation (96), the internal 
forces F can be expressed as a function of the displacement U. As an 
example, the values of F are represented for the 2-pixel of Figure 50 with 
A = l: 

Ft = C [(3 - 4i/)tt_i^ + (2 - 8i/)«a.i + (3 - 4u)ui,x + (10 - 8i/)«_i.o + (-36 + 48i/)«ojo 
+(10 - 8u)uija + (3 - 4i')u_i,-i + (2 - 8i/)«0, -1 + (3 - 4i/)«i,_i - 2v-i,i 
+2ui^ + 2i;_i,_i - 2«i,_i] 
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= C [(3 - 4i/)w_i^ + (10 - Bu)v^,i + (3 - 4j/)«i,i + (2 - 8i/)«_i,o + (-36 + 48i/)«o,o 
+(2 - 8i/)ui,o + (3 - 4i^)u_i,_a. + (10 - 8z/)uo.-i 4- (3 - 4i/)«i,_i - 2v-i^ 
+2«i,i + 2u_i,_i - 2t;i,_i] 



where : 



16(14- i/)(l-2«/) 



10 



and 



u 



-["1 



Equation (97) induces a linear relationship between a pixel and its 
neighbors. This relation is used to build the stiffness matrix of Equation (91). If 
(1=1 , 2) is considered as the displacement vector for the component Xj then: 



15 



where 



20 



iV* 



16(l+i/)(l-2i/) 



E 

16(l+i/)(l-2i/) 



3-4i/ 2-8i/ 3-4i/ 
10 -8z/ -36 + 48r/ 10 - 8i/ 
3-4z/ 2-8i/ 3-4i/ 



-2 0 2 
0 0 0 
2 0-2 
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16(l + i/)<l-2i/) 



-2 0 2 
0 0 0 
2 0-2 




16(1 -{- i^Xl - 2iy) 



3-4i/ 10-81/ 3-4z/" 

2- 8z/ -36+481/ 2-81/ 

3- 4i/ 10-81/ 3-4i/ 



The pairs CNl^,N'J , (i = 1, 2) will be referred to as the stiffness 
l<emels. 

Computation of the Displacement Vector 

The assumption made when calculating the. displacement vector will 
now be explained. Let v be a vertex of subdivided curve S and v' be Its 
corresponding vertex in the defonmed curve S' . Let us denote by U[v] the 
entry In the displacement vector U con-espondlng to v. Let us suppose that the 
displacement Is constant in each direction everywhere that is U[v] =.(ki, 
with ki, kae for ail v in S. Since the sum of all entries of either 
or Nljs zero, it follows that Fi[v] = F2M = 0 for all v In S which means that 
there is no internal force induced. The computation of the Intemai forces with 
the stiffness kernels is then invariant with respect to translation. 

Let vi, V2, V3, V4 and Vs be five adjacent vertices of S and v{ , v'^, Vj , , 
be their con-espondlng vertices in S' (see Figure 52). 

Let v,.,^be the X) coordinate of the spatial representation of the vertex Vi. 
Then, the following relation is obtained: 
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The translation invaiiance property leads to: 

where [v^^^ -Va^] stands for a matrix whose all entries equal v; ,_ -Vj^^ . The 

displacement component used to compute the internal force F| at vertex V3 is 
then: 

10 

which is the relative displacement of the vertex V3 with respect to V2.' However, 
nothing prevents computing the relative displacement of V3 with respect to V3. 
15 In this case, the Xk displacement component used to compute the internal 
force pi at vertex V3 would be: 

20 In order to take into account these facts, the average value of the 

relative displacements for all adjacent vertices to V2 is used. Thus: 

(97) 

25 
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Reorganizing the temis of Equation (97), for an arbitrary vertex V| of S, 
the following relation is obtained: 



L 2 J [ 2 



5 

or 



1 _ 8^S' '\ 

2 dtx^\ 



(98) 



If we assume a finite difference approximation of the second derivative of S 
1 0 and S' with respect to their vertices. 

Let us finally notice that the second derivative of the curve S in 
Equation (98) can be computed using Gaussian derivatives: 



15 dv^ 



where a controls the degree of smoothing. Such a computation of the second 
derivative of S allows to obtain smooth results by simulating a smoother target 
curve. 

20 

Experimental Results 



Active Contours 
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The approach proposed herein has been experimented on real and 
syntlietic Images in the context of high-resolution inriages of road databases. 
For each image, the results have been compared with another method. 

5 Figure 55a presents the results for the physics-based method (PB) 

according to the present invention for an aerial Image while Figure 55b shows 
the results for the finite element (FEM) method (a and p unknown). A material 
similar to rubber (see the Table In Figure 53) with E=150 and v=0.45 has 
been simulated. In both Images, the Image force Is the gradient (cr =1.5) of the 
10 bright line plausibility image obtained using a line detector proposed in [54] 
with the lirie detection scale set to 0.8. Figures 54a and 54b show respectively 
the initial curve So and the bright line plausibility image. 

Figure 56 shows a SAR image In which the initial curves are drawn in 
15 white. The line plausibility image (cr =1 .5 and line detection scale =0.8) of both 
bright and dark lines is shown in Figure 57. Figures 58 and 59 presents 
respectively the results obtained with PB (E=150 and v=0.45) and FEM (a 
and p unknown) methods. One can notice that the PB curves are closer to the 
shore than the FEM especially In region of high curvature. 
20 • 

Figure 60 shows some initialization for the first band of a Landsat 7 
image. The bright line plausibility image (<t=1.5 and line detection scale =0.8) 
is shown in Figure 61. Figures 62 and 63 present the results obtained with the 
PB and FEM methods respectively. 

25 

The fact that the deformations obtained using the model according to 
this illustrative embodiment of the present Invention have a physical 
interpretation has been discussed. The fact that the objects modeled using the 
PB method have their own physical properties and the ability to recover their 
30 original shape when the external forces applied on them are removed has also 
been discussed. To Illustrate this fact, Figures 64a and 64b show some 
Initialization and the corrected curve for a synthetic Image. Figures 65a 
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through Figures 65d show the evolution of this corrected curve when the 
external forces are removed. Figure 65d presents both the final curve (in 
black) and the Initial curve (in white). One can clearly notice that the curve has 
recovered its original shape but has also experienced a spatial shift. 

5 

Conclusion 

A new model for active contours was presented. The proposed 
approach decompose the Image using an image model based on algebraic 
10 topology. This model uses generic mathematical tools which can be applied to 
solve other problems such as linear and non-linear diffusion and optical flow 
[57]. Moreover, the model wori<s with either 2D or 3D Images and can easily 
be extended to active surfaces and active volumes. 

15 The approach presents the following differences with the other 

nnethods: 1) Both global and local quantities are used; 2) The model is based 
upon basic laws of physics. This allows us to give a physical explanation to the 
deformation steps; 3) The curves and surfaces have physical behaviors such 
as the ability to recover their original shape once the applied forces are 

20 removed; 4) Approximation are made only when the constitutive equation is 
involved. 

Although the present invention has been described hereinabove by way 
of non-restrictive illustrative embodiments thereof, it can be modified at will, 
25 within the scope of the appended claims, without departing from the spirit and 
nature of the subject invention. 
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WHAT IS CLAIMED IS: 

1. A cx)mputational image model, comprising: 

an Image support including a structure of n-pixels comprising pixel 

faces; 

quantities related to image features; and 

an algebraic stmcture relating the quantities to the n-plxels and/or pixel 
faces, the algebraic structure comprising algebraic operations defining a 
relation between the quantities.% 

2, A computational image model as defined In claim 1, wherein each n- 
pixel is defined as a geometrical structure comprising vert:lces, edges, faces 
and a volume, and wherein each n-pixel comprises: 

a first pixel dimension n=0 including the vertices of the n-pixel; 
a second pixel dimension n=1 including the edges of the n-pixel; 
a third pixel dimension n=2 including the faces of the n-pixel; 
a fourth pixel dimension n=3 including the volume of the n-pixel; and 
a n* pixel dimension n including the hypen/olume of the n-pixel. 

3, A computational image model as defined In claim 1, wherein the 
geometrical structure is selected from the group consisting of: a cube, a 
triangle, a hexagone and a pentagone. 

4. A computational image model as defined in claim 1, wherein the 
quantities related to image features are selected from the group consisting of: 
scalar quantities, vectors, tensors and matrices. 

5, A computational image model as defined in claim 1, wherein the 
algebraic operations comprise problem-independent operations. 

6. A computational image model as defined in claim 1, wherein the 
algebraic operations comprise problem-dependent operations. 
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7. A computational image model as defined in claim 1, wherein the 
structure of n-pixels comprises pairs of disjoint n-pixels. 

5 8. A computational image model as defined In claim 1, wherein the 

structure of n-pixels comprises pairs of n-pixels intersecting through a common 
i-pixel, where i < n. 

9. A computational image model as defined in claim 1 . wherein each n- 
10 pixel is translated algebraically into a q-pixel, wherein q e {1, 2,.,., n}. 

10. A computational image model as defined in claim 9, wherein each 
q-pixel includes (q-1)-faces, (q-2)-faces. .... (q-q)-faces. 

15 11. A computational image model as defined In claim 9, wherein the 

image support comprises a geometrical complex, which is a collection of q- 
pixels. 

12. A computational image model as defined in claim 10, wherein the 
20 image support comprises a geometrical complex, which is a collection of q- 

pixels, and wherein: 

- every face of a q-pixel in the geometrical complex is also located in the 
geometrical complex; and 

- any pair of two q-pixels of the geometrical complex have an intersection 
25 which is either empty or constituted by a common face of both q-pixels 

of the pair. 

13. A computational image model as defined in claim 11, comprising a 
plurality of image supports forming the geometrical complex. 

30 
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14. A computational image model as defined in claim 11, wherein the 
geometrical complex is expressed in algebraic form as a q-chain, which Is a 
linear combination of all the q-pixels of the geometrical complex. 

15. A computational image model as defined in claim 9, wherein the 
geometrical complex comprises q-cochains, which are relations associating 
quantities related to image features to the q-pixels and/or faces of said q- 
pixels. 

16. A computational Image model as defined In daim 15. wherein the 
quantities related to image features and associated to the q-plxels and/or 
faces of said q-pixels are global quantities associated to all the q-pixels. 

17. A computational image model as defined in claim 15, wherein the 
quantities related to image features and associated to the q-plxels and/or 
faces of said q-pixels are local quantities each associated to one q-pixel and/or 
faces of said one q-pixei. 

18. A computational image model as defined in claim 16, comprising 
(q>1 )-cochains to represent the local quantities. 

19. A computational image model as defined in claim 17, comprising 0- 
cochain to represent the global quantities. 

20. A computational image model as defined in claim 17, wherein the 
algebraic operations comprise a coboundary operation giving a relationship 
between the q-cochains. 

21. A computational image model as defined in claim 9, wherein: 

the image support comprises a plurality of geometrical complexes, 
each being a collection of q-pixeis; and 
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the algebraic operations comprise a codual operation establishing a link 
between q-cochains that belong to different geometrical complexes. 

22. A method of computationally modelling an image, comprising: 
producing an Image support including a structure of n-pixels comprising 

pixel faces; 

defining quantities related to image features; and 

relating the quantities to the n-pixels and/or pixel faces through an 
algebraic structure, and relating the quantities to each other through algebraic 
operations. 

23. A method of computationally modelling an Image as defined in 
claim 22, wherein relating the quantities to the n-pixels and/or pixel faces 
through an algebraic structure comprises translating each n-plxel algebraically 

Into a q-pixel, wherein q e {1,2 n}, wherein each q-pixel includes (q-1)- 

faces, (q-2)-faces (q-q)-faces. 

24. A method of computationally modelling ^n image as defined in 
claim 22, wherein producing an image support comprises fomiing a 
geometrical complex, which Is a collection of q-pixels, and wherein: 

- every face of a q-pixel In the geometrical complex is also located In the 
geometrical complex; and 

- any pair of two q-pixels of the geometrical complex have an intersection 
which Is either empty or constituted by a common face of both q-pixels 
of the pair. 

25. A method of computationally modelling an image as defined in 
claim 24, wherein producing an image support comprises forming a plurality of 
image supports fomning the geometrical complex. 

26. A method of computationally modelling an Image as defined in 
claim 24, wherein relating the quantities to the n-pixels and/or pixel faces 
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through an algebraic structure comprises expressing the geometrical complex 
in algebraic form as a q-chain, which is a linear combination of all the q-pixels 
of the geometrical complex. 

27. A method of computationally modelling an image as defined in 
claim 24, wherein relating the quantities to the n-pixels and/or pixel faces 
through an algebraic structure comprises forming, in the geometrical complex, 
q-cochains which are relations associating quantities related to image features 
to the q-pixels and/or faces of said q-pixels. 

28. A method of computationally modelling an image as defined in 
claim 22, wherein defining quantities related to Image features comprises 
defining global quantities associated to all the q-pixels. 

15 29. A method of computationally modelling an image as defined In 

claim 22, wherein defining quantities related to image features comprises 
defining local quantities associated to one q-pixel and/or faces of said one q- 
pixel. 

20 30. A method of computationally modelling an image as defined in 

claim 27, wherein relating the quantities to each other . through algebraic 
operations comprise producing a coboundary operator giving a relationship 
between q-cochains. 

25 31. A method of computationally modelling an image as defined in 

claim 27, wherein: 

producing an image support comprises forming a plurality of 
geometrical complexes, each being a collection of q-pixels; and 

relating the quantities to each other through algebraic operations 
30 comprises producing a codual operation establishing a link between cochains 
that belong to different geometrical complexes. 
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32. An image modelling method as defined in claim 27, wherein relating 
the quantities to the n-pixels and/or pixel faces through an algebraic structure 
comprises expressing a global quantity associated with all q-pixels through a 
q-cochain such that, for two adjacent q-plxels and c^, the q-cochain 
satisfies the relation i\c\ + \cl) = \F^ (c\) + X,F^ {cD , where X^ and X2 are 
integers. 

33. An image modelling method as defined in claim 22, wherein: 

- relating the quantities to the n-pixels and/or pixel faces through an 
algebraic structure comprises translating each n-pixel algebraically into a 
q-pixel, wherein q e {1, 2,..., n}, wherein each q-pixel includes (q-l)-faces, 
(q-2Haces, {q-q)-faces; 

- producing an image support comprises forming geometrical complexes, 
each being a collection of q-plxels; 

- relating the quantities to the n-plxels and/or pixel faces through an 
algebraic structure comprises: 

o expressing each geometrical complex in algebraic fonn as a q- 

chain, which is a linear combination of all the q-plxels of the 

geometrical complex; 
o forming, in the geometrical complexes, q-cochains which are 

relations associating quantities related to image features to the 

q-pixels and/or faces of said q-pixels; 

- relating the quantities to each other through algebraic operations 
comprises: 

o producing a coboundary operator giving a relationship between 

the q-cochains; and 
o producing a codual operation establishing a link between q- 

cochalns that belong to different geometrical complexes. 

34. A cx)mputational framewori< for solving a problem using an image 
computationally modelled by means of the method of claim 33, comprising: 
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identifying basic laws associated to the problem; 
from the identified basic laws, defining quantities related to the problem; 
associating the quantities to respective q-cochains; 
associating the basic laws related to the problem to respective 
5 coboundary and codual operations; and 

resolving the resulting algebraic system. 

35- A computational framework as defined In claim 34, wherein forming 
geometrical corhplexes comprises forming first and second geometrical 
10 complexes. 

36. A computational framework as defined in claim 35. wherein 
Identifying basic laws associated to the problem comprises supporting one 
basic law through the first geometrical complex. 

15 

37. A computational framework as defined In claim 36, wherein the 
problem to be solved Is a 2D global differential equation for heat flow in a 
homogeneous medium, and wherein said one basic law is a heat flow law. 

20 38. A computational framewori< as defined in claim 37, wherein 

associating the quantities to respective q-cochains comprises representing a 
global quantity of temperature through a 0-cochain, and associating the heat 
flow law through a 1-cochain. 

25 39. A computational framework as defined in claim 35, wherein 

identifying basic laws associated to the problem comprises supporting one 
basic law through the second geometrical complex. 

40. A computational framework as defined in claim 39, wherein the 
30 problem to be solved is a 2D global differential equation for heat flow in a 
homogeneous medium, and wherein said one basic law is a heat source law. 
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41. A computational framework as defined in claim 36, wherein 
identifying basic laws associated to the problem comprises supporting a 
second basic law through the second geometrical cornplex, and wherein 
associating the basic laws related to the problem to respective coboundary 

5 and codual operations comprises representing a constitutive law linking basic 
laws from the first and second geometrical complexes by a codual operation. 

42. An image modelling method as defined in claim 22, wherein: 
relating the quantities to the n-pixels and/or pixel faces through an 
algebraic structure comprises translating each n-pixel algebraically into a 
q-pixel, wherein q e {1, 2...., n}, wherein each q-pixel includes (q-l)-faces, 
(q-2)-faces, (q-q)-faces; 

producing an image support comprises forming a geometrical complex, 
which is a collection of q-pixels; 

relating the quantities to the n-pixels and/or pixel faces through an 
algebraic structure comprises: 

o expressing the geometrical complex in algebraic form as a q- 
chain, which is a linear combination of all the q-pixels of the 
geometrical complex; 
o forming, in the geometrical complex, q-cochains which are 
relations associating quantities related to image features to the 
q-pixels and/or faces of S9id q-pixels; 
relating the quantities to each other through algebraic operations 
comprises: 

o producing coboundary operations giving a relationship between 
the q-cochains. 

43. A computational framework for solving a problem using an image 
computationally modelled by means of the method of claim 42, comprising: 

30 identifying basic laws associated to the problem; 

from the identified basic laws, defining quantities related to the problem; 
associating the quantities to respective q-cochains; 
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associating the basic laws related to the problem to respective 
coboundary operations; and 

resolving the resulting algebraic system. 

44. A computational framework for solving a heat transfer problem, 
comprising: 

producing an image support including a structure of n-pixels, the image 
support comprising: 

6 q-pixels respectively translating the n-pixel algebraically, 

wherein q e {1, 2 n}, and wherein each q-pixel includes (q- 

1)-faces, (q-2)-faces, (q-q)-faces; 
o geometrical complexes each being a collection of q-plxels; 
o q-chains respectively expressing the geometrical complexes in 
algebraic form, each q-chain being a linear combination of all 
the q-pixels of the geometrical complex; 
o in the geometrical complexes, q-cochains which are relations 
associating quantities related to image features to tine q-pixels 
and/or faces of said q-pixels; and 
o a coboundary defining a relation between q-cochains; 
computing a q-coghain T of a first of said geometrical complexes as the 
location of unl<nown temperatures; 

computing a q-cochain H of the first geometrical complex as a global 
temperature variation; 

finding a q-cochain s of a second geometrical complex as a global 
energy variation, as a function of the q-cochain H through a linear 
transformation; 

finding the q-cochain 8 as a function of the q-cochain T; 
defining a q-cochain G of the first geometrical complex from the q- 
cochain T through a first coboundary operation, transforming the q-cx)chain G 
into a q-cochain Q of the second geometrical complex, and defining, from the 
q-cochain Q and through a second coboundary operation, a q-cochain D of the 
second geometrical complex as a global diffusion; 
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defining a q-cochain S of the second geometrical complex as a global 
source; and 

establishing a relation between the q-cochains s, D and S. 

45, A computational framework for two-dimensional active contour 
model, comprising: 

producing an image support including a structure of n-pixels, the image 
support comprising: 

- q-pixels respectively translating the n-pixel algebraically, 
wherein q e {1, 2,..., h}. and wherein each q-pixel includes 
(q-l)-faces, (q-2Haces, (q-q)-faces; 

- geometrical complexes each being a collection of q-pixels; 

- q-chains respectively expressing the geometrical complexes 
in algebraic form, each q-chain being a linear combination of 
all the q-pixels of the geometrical complex; 

- in the geometrical complexes, q-cochains which are relations 
associating quantities related to image features to the q- 
pixels and/or faces of said q-pixels; and 

- a coboundary defining a relation between q-cochains; 
computing a displacement q-cochain D of a first of said geometrical 

complexes; 

computing a strain q-cochain S . of a second of said geometrical 
complexes, comprising: 

- defining an approximate strain function s{x) as a function of 
the q-cochain D; 

- expressing the q-cochain S as a function of the approximate 
strain function and relative positions of the first and second 
geometrical complexes; and 

computing a force q-cochain F of the second geometrical complex as a 
coboundary of the strain q-cochain S. 
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SUPPORT ALGEBRAICALLY IN RELATION TO 
THE n-PIXEL DIMENSIONS 
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PRODUCING A GEOMETRICAL COMPLEX 
INCORPORATING AT LEAST ONE n-PIXEL 
IMAGE SUPPORT 






EXPRESSING THE GEOMETRICAL COMPLEX 
ALGEBRAICALLY BY MEANS OF q-CHAINS 






EXPRESSING THE SCALAR VECTORS AND 
MATRICES BY MEANS OF THE COEFFICIENTS 
OF THE q-CHAIN 






EXPRESSING GLOBAL QUANTITIES WITH ALL 
q-PIXELS THROUGH q-COCHAINS Fq 
COMPRISING ASSOCIATING GLOBAL 
QUANTITIES TO THE q-PIXELS AND ITS FACES 
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QUANTITIES RELATED TO THE PROBLEM 




r 



6701 



6702 



6703 



THE DEFINED QUANTITIES AND 
IMAGE SUPPORT THROUGH AN ALGEBRAIC 
LANGUANGE. INCLUDING INSTANTIATION OF 
THE COBOUNDARY AND CODUAL 
OPERATORS 



6704 



1 


r 


RESOLUTION OF THE RESULTING ALGEBRAIC 
SYSTEM 



6705 



il 



(12) DMTERNAITONAL APPLICATION PUBLISHED UNDER THE PATENT COOPERATION TREATY (PCT) 

ill 



(19) World InteUectual Property 
Organization 
Internationa] Bureau 




(43) International Publication Date 
19 February 2004 (19.02.2004) PCT 



(10) International Publication Number 

WO 2004/015631 A3 



(51) International Patent Classification^: G06T 5/00, 

9/00. G06F 17/13 

(21) International Application Number: 

PCT/CA2003/001214 

(22) International Filing Date: 8 August 2003 (08.08.2003) 

(25) Filing Language: English 

(26) Publication Language: English 



(30) Priority Data: 

2,397,389 



9 August 2002 (09.08.2002) CA 



(71) Applicants (for all designated States except US): UNI- 
VERSITE DE SHERBROOKE [CA/CA]; 2500, boule 
vard de ITJniversite, Sherbrooke, Quebec JIK 2R1 (CA). 
BISHOP'S UNIVERSITY [CAJCA]; LennoxviUe, Que- 
bec JIM 1Z7 (CA). 



(72) Inventors; and 

(75) Inventors/AppOcants (for US only): ZIOU, DJemel 
[CA/CA]; 2665 Maricourt, Sherbrooke. Quebec JIK 1R6 
(CA). AUCLAIR-FORTIER, Marie-Flavie [CA/CA]; 
264 Lockwell, Quebec, Quebec GIR 1 V9 (CA). ALLILI, 
IVladJid [CA/CA]; 585 London, apt. 1, Sherbrooke. 
Quebec JIH 3N2 (CA). 

(74) Agents: BROUILLETTE, Robert et al.; Brouillette 
Kosie Prince, 1100 Rene-Levesque Blvd. West, 25th 
Floor, Montreal. Quebec, H3B 5C9 (CA). 

(81) Designated States (national): AE, AG, AL, AM, AT, AU, 
AZ, BA, BB, BG, BR, BY, BZ, CA, CH, CN, CO, CR, CU, 
CZ, DE, DK, DM, DZ, EC, EE, ES, H, GB, GD, GE, GH, 
GM, HR, HU, ID. IL, IN, IS. JP. KE, KG. KP. KR. KZ, LC, 
LK, LR, LS, LT. LU, LV, MA. MD, MG, MK, MN, MW, 
MX, MZ. NT. NO. NZ. OM, PG, PH, PL, FT, RO, RU, SC, 
SD, SE, SG, SK, SL. SY, TJ. TM. TN, TR, TT. TZ, UA, 
UG, US. UZ. VC, VN, YU. ZA. ZM, ZW. 



f Continued on next page] 



(54) Title: IMAGE MODEL BASED ON N-PDCELS AND DEFINED IN ALGEBRAIC TOPOLOGY AND APPLICATIONS 
THEREOF 



DEFINING A PIXEL DIMENSION n TO PRODUCE 
A n-PIXEL IMAGE SUPPORT 




r 


EXPRESSING THE n-PIXELS OF THE IMAGE 
SUPPORT ALGEBRAICALLY IN RELATION TO 
THE n-PIXEL DIMENSIONS 







6601 



6602 



PRODUCING A GEOMETRICAL COMPLEX 
INCORPORATING AT LEAST ONE n-PIXEL 
IMAGE SUPPORT 



(57) Abstract: A computational image model comprises an image 
support including a structure of n-pixels comprising pixel faces, 
quantities related to image features, and an algebraic structure 
relating the quantities to the n-pixels and/or pixel faces, the algebraic 
structure comprising algebraic operations defining a relation 
between the quantities. A method of computationally modelling an 
image comprises producing an image support including a structure 
of n-pixels comprising pixel faces, defining quantities related to 
image features, and relating the quantities to the n-pixels and/or 
pixel faces through an algebraic structure, and relating the quantities 
to each other through algebraic operations. 



6603 



ID 



o 



EXPRESSING THE GE 
ALGEBRAICALLY BY Ml 


1 r 

©METRICAL COMPLEX 
EANS OF q-CHAINS 






EXPRESSING THE SCALAR VECTORS AND 
MATRICES BY MEANS OF THE COEFFICIENTS 
OF THEq-CHAiN 






EXPRESSING GLOBAL QUANTITIES WITH ALL 
q-PlXELS THROUGH q-COCHAINS Fq 
COMPRISING ASSOCIATING GLOBAL 
QUANTTTIES TO THE q-PIXELS AND ITS FACES 



6604 



6605 



6606 



0 2004/015631 A3 liliilililiillilllliliiililliiiiliH 



(84) Designated States (regional): ARIPO patent (GH, GM, 
KE, LS, MW, MZ, SD, SL, SZ, TZ, UG, ZM, ZW), 
Eurasian patent (AM, AZ, BY, KG, KZ, MD, RU, TJ, TM), 
European patent (AT, BE, BG, CH, CY, CZ, DE, DK, EE, 
ES. H. FR, GB. GR. HU. IE. IT. LU. MC. NL, PT, RO. 
SE, SI, SK, TR), OAPI patent (BF, BJ, CF, CG, CI, CM, 
GA, GN, GQ, GW. ML, MR. NE, SN, TD. TG). 

PubUshed: 

— with international search report 



— before the expiration of the time limit for amending the 
claims and to be republished in the event of receipt of 
amendments 

(88) Date of publication of tlie international search report: 

22 July 2004 

For two 'letter codes and other abbreviations, refer to the "Guid- 
ance Notes on Codes and Abbre\nations" appearing at the begin- 
ning of each regular issue of the FCT Gazette, 



INTiH^ATIONAL SEARCH REPORT 



A. CLASSIFICA1IONOF. 

IPC 7 G06T5/00 



MATTER , 

G06T9/00 




l-AppOcation No 

'PCT/CA 03/01214 



606F17/13 



According to International Patent ClassificaKon (IPC) or to both natlonai dasslBcalfon and IPC 



B. RELDS SEARCHED 


Minimum documentation searched (das^fication system followed by classification symbols) 

IPC 7 606T G06F 


Documentation searched other than nninimum documentation to the extent that such documents are included In the fields seaiched 


Electronic data base consulted during the International search (name of data bas 


e and, where practical, search terms used) 


EPO-Internal , WPI Data, PAJ, INSPEC 






C. DOCUMENTS CONSIDERED TO BE RELEVANT 


Category • 


Citation of document, with indication, where appropriate, of the relevant passages 


Relevant to daim No. 


X 
A 


ZIOU D, ALLILI M: "An Image model vilth 
roots In computational algebraic topology: 
A primer" 

TECHNICAL REPORT 264, 

April 2001 (2001-04), XP001189245 
Departei|ient de mathematlques et 
d'informatlque, University de Sherbrooke, 
SherbrooKe, Canada 

* abstract * 

* sections 2, 4, 5, 6.1 * 


1-43 
44,45 






/— 




1 X| Further documents are listed in the continuation of box C. 


1 j Patent family memtiers are listed In annex. 


Special categories of cited documents : 

"A" document defining the genera] state of the art which is not 
considered to be of particular relevance 

'E* earlier document but published on or after the international 
filing dale 

*L* document which may throw doubts on priority ctalm(s) or 
which Is cited to establish the publication date of another 
citation or other special reeison (as specffied) 

"O" document referring to an oral disclosure, use, exhlbftlonor 
other means 

'P' docunrteht published prior to the International fifing date but 
later than the priority date claimed 


"T" later document published after the international filing date 
.or priority date and not in conflict with the application but- 
cited to understand the principle or tiieory underiylng the 
Invention 

'X' document of particular relevance; the claimed invention 
cannot be considered novel or cannot be considered to 
involve an Inventive step when the document is taken alone 

'Y* document of particular relevance; the dafmed Invention 

cannot be considered to involve an Inventive step when the 
document is combined with one or more otiier such docu- 
ments, such combination being obvious to a person sidlled 
in the art 

'&* document memberofthe same patent family 


Data of the 


actual oompletton of the international search 


Date of mailing of the Intemallona! sea 


liGh report 


12 May 2004 


26/05/2004 




Name and nnalfing address of the ISA 

European Patent Office, P.a 5818 Palentiaan 2 
NL-2280HVRqswl3lc 
TeL (+31-70) 340-^)40. Tx. 31 651 epo nl, 
Fax (431-70) 340-3016 


Authorized officer 

. Domingo Vecchloni 





Form PCT/1SA/210 (second sheat) (January 2004) 



# 



intiBational search report 



Int^H^onal Application No 

PCT/CA 03/01214 



C.(CoiitlnuaUon) DOCUMENTS CONSIDEREO TO BE RELEVANT 



Category " 



CBalton of document, with Indication, wtiere appropriate, of the relevant passages 



Relevant to cteim No. 



AUCLAIR-FORTIER M F, POULIN P, ZIOU D, 
ALLILI M: "Physics-based resolution of 
diffusion and optical flow: A 
computational algebraic topology approach" 
TECHNICAL REPORT 269, 

November 2001 (2001-11), XP001189246 
Departement de mathematiques et 
d'informatique, Universlte de Sherbrooke, 
Sherbrooke, Canada 
cited in the application 

* abstract * 

* sections 4, 5 * 

POULIN P, AUCLAIR-FORTIER M F, ZIOU D, 
ALLILI M: "A physics-based model for the 
deformation of curves: A computational 
algebraic topology approach" 
TECHNICAL REPORT 270, 
8 February 2002 (2002-02-08), XP001189247 
Departement de mathematiques et 
d'informatique. University de Sherbrooke, 
Sherbrooke, Canada 
cited in the application 

* abstract * 

* sections 3-5 * 

POULIN P, AUCLAIR-FORTIER M F, ZIOU D, 
ALLILI M: "A physics-based model for 
active contours: A computational algebraic 
topology approach" 

SYMPOSIUM ON GEOSPATIAL THEORY, PROCESSING 
AND APPLICATIONS, OTTAWA, CANADA, 
^Online! 9-12 July 2002, XP002279512 
Retrieved from the Internet: 
<URL : http : //www, 1 sprs . org/commi ssi on4/proc 
eedi ngs/pdf paper s/375 . pdf > . „ 

* retrieved on 2004-05-11! 

* section 3, paragraph 1; references * 



1-45 



1-45 



Foan PCTyiSA/210 (oonUnuatlon of second ^leet) (Januaiy aXM) 



This Page is Inserted by IFW Indexing and Scanning 
Operations and is not part of the Official Record 

BEST AVAILABLE IMAGES 

Defective images within this document are accurate representations of the original 
documents submitted by the applicant. 

Defects in the images include but are not limited to the items checked: 

□ BLACK BORDERS 

□ IMAGE CUT OFF AT TOP, BOTTOM OR SIDES 
-□TfaDED TEXT OR DRAWING 

□ BLURRED OR ILLEGIBLE TEXT OR DRAWING 

□ SKEWED/SLANTED IMAGES 

□ COLOR OR BLACK AND WHITE PHOTOGRAPHS 

□ GRAY SCALE DOCUMENTS 

-''Q^INES OR MARKS ON ORIGINAL DOCUMENT 

□ REFERENCE(S) OR EXHIBIT(S) SUBMITTED ARE POOR QUALITY 

□ OTHER: 

IMAGES ARE BEST AVAILABLE COPY. 
As rescanning these documents will not correct the image 
problems checked, please do not report these problems to 
the IFW Image Problem Mailbox. 



